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ABSTRACT 


A  numerical  model  to  simulate  the  deformations  in  gear  teeth  subjected  to 
rolling  loads  as  in  the  ausrolling  process  has  been  developed.  Ausrolling  involves 
applying  rolling  loads  to  the  gear  when  it  is  in  the  metastable  austenitic  state.  A 
model  of  the  process  must  consider  material,  geometric  and  surface  non-linearities 
as  well  as  changes  in  temperature  and  material  properties  with  time  and  rolling 
loads  in  three  dimensions.  Only  some  of  these  requirements,  namely,  the  elastic- 
plastic  flow;  geometric  non-linearities  due  to  large  deformations;  frictional  contact 
conditions  at  the  die-workpiece  interface;  and  the  travelling  loads  due  to  rolling, 
have  been  considered  here  to  be  of  primary  importance.  The  objective  of  this  thesis 
is,  accordingly,  to  satisfactorily  establish  the  foremen tioned  features  in  the  non¬ 
linear  finite  element  analysis  program  (NOFEAP).  The  theoretical  aspects  of  the 
non-linear  formulations  have  been  briefly  described  and  the  model  of  the  rolling 
process  has  been  outlined  here.  The  implemented  non-linear  formulations  have 
been  compared,  individually  and  in  various  combinations,  with  existing  analytical, 
numerical  and  experimental  solutions.  Since  no  such  results  are  available  in  the 
literature  to  compare  the  response  of  the  rolling  model,  several  experiments  have 
been  conducted  on  aluminum  6061  disks.  The  experimental  results  and  their 
numerical  counterparts  have  been  included  with  discussion  on  the  similarities  and 
differences  in  the  two.  Finally,  the  applicability  of  NOFEAP  to  modelling  the 
deformations  in  gear  rolling  has  been  demonstrated  by  modelling  the  deformations 
in  the  approach  and  the  trailing  sides  of  a  gear  tooth  when  rolled  against  a  hard 
die.  These  deformation  plots  and  stress  contours  are  also  given  in  this  thesis. 
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Chapter  1 
INTRODUCTION 


A  low  temperature  thermomechanical  process  is  being  developed  by  the  author 
and  coworkers  for  finishing  precision  machine  elements.  This  process  is  being 
applied  to  the  finishing  of  gears  by  working  the  surface  layers  of  the  teeth  by  a 
combination  of  rolling,  swaging  and  shearing  operations.  This  process  is  referred 
to  as  “ausrolling”  and  results  in  a  significant  saving  in  time  and  processing  costs 
over  conventional  finishing  operation  such  as  grinding.  Ausrolling  also  results  in  a 
fine  grained  martensitic  structure  and  high  compressive  stresses  in  the  carburized 
surface  which  improves  the  fatigue  life  of  the  machine  elements.  One  key  to  the 
success  of  ausrolling  is  the  ability  to  reliably  and  accurately  predict  the  finished 
dimensions  of  the  machine  components.  Such  an  estimate  can  be  obtained  by 
numerically  analyzing  the  process.  This  thesis  presents  a  non-linear  finite  element 
model  used  to  study  the  ausrolling  process  which  is  described  next. 

1.1  The  Ausrolling  Process 

Most  mechanical  systems  contain  a  few  precision  machine  elements  which  must 
possess  superior  mechanical  properties  to  improve  performance,  and  they  must  be 
precisely  finished  to  minimize  noise  and  vibration.  These  requirements  are  satisfied 
by  a  number  of  manufacturing  steps,  as  shown  in  Figure  1.1.  Unfortunately, 
some  of  these  steps  may  annul  one  or  more  of  the  desirable  characteristics 
introduced  by  the  preceding  operations.  Moreover,  as  the  specifications  for  the 
mechanical  properties  and  dimensional  tolerances  become  more  stringent,  the 
cost  of  processing  can  be  significantly  higher.  These  detrimental  aspects  of  the 
conventional  finishing  processes  can  be  neutralized  if  some  of  the  operations  can  be 
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merged  into  a  single  process.  Such  a  combination  of  steps  has  been  accomplished 
in  the  ausrolling  process  (Figure  l.l). 

Since  the  plastic  deformations  are  introduced  by  ausrolling  in  the  machine 
component  while  the  component  is  in  metastable  austinitic  state,  the  time- 
temperature-transformation  (TTT)  diagram  of  the  material  plays  an  important 
role  in  the  operation.  There  are  two  essential  requirements  of  a  steel  to  be  effective 
in  an  ausforming  operation:  presence  of  some  carbide  forming  elements  and  a 
deep  austenite  bay  region  in  the  TTT  curve.  The  carburized  9310  steel  (1.0C- 
1.2Cr-3.25Ni-0.13Mo),  which  has  been  successfully  surface- ausformed,  follows  the 
TTT-diagram  shown  in  Figure  1.2.  The  mechanical  system  used  for  ausforming  is 
described  next. 

Ausrolling,  applied  to  gears  [lj,  is  a  modified  form  of  the  conventional  gear 
rolling  operation,  where  the  die  gear  (A)  is  powered  and  drives  the  softer 
workpiece  gear  (B),  as  shown  in  Figure  1.3.  The  workpiece,  henceforth  referred 
to  as  the  gear,  is  fed  into  the  die  by  a  combination  of  horizontal  and  vertical 
movements.  This  operation,  where  the  gear  is  simultaneously  worked  along  the 
circumferential  surface  and  the  thickness,  is  called  the  swage-rolling  process.  The 
slide  (C)  on  which  the  gear  is  mounted  and  the  rear  column  (D)  are  connected  to 
electrohydraulically  operated  actuators.  The  vertical  feed  is  achieved  by  lowering 
the  front  ram,  thus  moving  the  gear  across  the  die.  The  horizv..cai  feed  is  obtained 
by  translating  the  vertical  motion  of  the  rear  ram  into  a  horizontal  motion  by 
the  use  of  a  40:1  tapered  column  (D)  between  the  rear  plate  (E)  and  a  block  (F) 
wedged  against  the  middle  plate  (G).  An  indexing  gear  (H)  may  be  used  to  enforce 
pure  rolling  at  desired  radius  and  to  assure  smooth  entry  into  the  die.  The  whole 
assembly  is  immersed  into  a  hot  oil  tank  to  maintain  temperature  and  lubrication 
conditions.  The  loading  conditions  range  from  pure  rolling  at  the  pitch-line  to 
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combined  rolling  and  sliding  away  from  the  pitch-line.  Thus  the  state  of  stress 
changes  radically  along  the  active  profile  of  the  involute  gear  surface.  Material  on 
the  approach  side  of  the  workpiece  gear  is  displaced  toward  the  pitch-line  while  it 
is  displaced  in  the  opposite  direction  on  the  trailing  side  as  shown  in  Figure  1.4. 

GEAR  R OTATI  ON 


(driven) 


DIE  ROTATION 
(driver) 

Ficur-  1.4  RELATIVE  SLIDING  BETWEEN  GEAR  TEETH  IN  MESH 

In  addition  to  stress  gradients  on  the  surface  of  the  tooth  being  worked  there  are 
also  compositional  and  thermal  gradients  in  the  direction  normal  to  the  surface. 
The  compositional  gradients  are  due  to  the  carburized  condition  wniie  the  thermal 
gradients  are  a  result  of  the  quenching  to  the  metastable  austenitic  state  just 
prior  to  and  during  working.  Significant  dynamic-hardening  is  also  characteristic 
of  the  process  since  working  takes  piace  beiow  the  recr/stailization  temperature 
and  under  conditions  that  promote  rapid  dislocation  network  pinning.  When  this 
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process  is  applied  to  spur  gears  a  specially  tapered  lead  geometry  is  provided  In 
the  die  (Figure  1.5)  to  produce  a  swaging  action  on  the  teeth  during  forming.  The 
lead  in  the  die  being  used  in  the  experiments  also  contains  a  hollow  section  to 
produce  a  crown  contour  in  the  gear. 


GEAR 

ENTRY 


Figure 


1.5  LEAD  PROFILE  FOR  GEAR  DIE 


1.2  Scope  and  Outline  of  the  Thesis 

A  numerical  model  which  adequately  simulates  forming  processes  must  account 
for  material,  geometric  and  surface  non-linearities  depending  on  the  application. 
Formulations  of  these  non-linearities  abound  in  the  literature  [2—4}.  The  objective 
of  the  thesis  is  to  combine  these  non-linearities  in  one  algorithm  so  that  the 
resulting  model  can  be  used  to  study  the  ausrolling  process. 

Important  steps  of  the  non-linear  formulations  are  described  in  Chapter  2. 
Each  formulation  is  preceded  by  a  review  of  the  available  literature  on  the 
subject.  Material  non-linearities  are  modelled  by  an  elastic-plastic  formulation 
[5j.  Geometric  non-linearities  are  formulated  by  using  the  total  and  the  updated 
Lagrangian  formulations  [3j.  The  surface  non-linearities  are  incorporated  into 
the  model  by  using  a  frictional  contact  formulation  [6].  The  transformation 
of  these  formulations  into  the  finite  element  matrix  equations  [7]  and  a  brief 
description  of  the  non-linear  finite  element  analysis  program  (NOFEAP)  [8]  which 
combines  these  models  into  a  computer  code  are  also  presented  in  Chapter  2. 
The  implementation  of  each  of  these  non-linear  formulations  has  been  verified  by 
solving  a  number  of  benchmark  problems  using  NOFEAP.  The  results  of  these 
verification  excercises  are  presented  in  Chapter  3. 

The  most  unique  feature  of  the  ausrolling  process  is  that  the  deformation  in 
the  workpiece  is  introduced  by  rolling  the  gear  against  a  hard  die.  The  rolling 
process  has  not  been  widely  modelled  in  existing  finite  element  codes.  Mori  [9] 
has  used  rigid-plastic  finite  element  method  to  study  the  rolling  of  sheet  metal 
by  approximating  the  tractions  at  the  roll-strip  interface  and  also  the  velocity  of 
the  roll  relative  to  the  strip  at  its  neutral  region.  Elastic-plastic  rolling  contact 
has  been  studied  by  Bhargava  [10,11]  who  simulates  the  rolling  conditions  by 


translating  a  semi-elliptic  Hertzian  pressure  distribution  along  the  surface  of 
a  flat  workpiece.  These  methods  which  assume  a  pressure  distribution  have 
severe  drawbacks  when  applied  to  rolling  cylinders  or  gears.  A  model  of  the 
rolling  process  developed  using  geometric  and  kinematic  considerations  is  used  in 
NOFEAP.  The  details  of  this  model  and  its  application  to  the  study  of  rolling 
cylinders  and  gear  teeth  are  presented  in  Chapter  4.  The  conclusions  of  this  thesis 
are  discussed  in  Chapter  5. 

1.3  Review  of  Literature 

There  is  avast  amount  of  literature  on  non-linear  finite  element  analysis.  For  the 
sake  of  organized  presentation  the  literature  review  is  divided  into  various  sections. 
This  section  only  summarizes  the  research  on  combined  non-linear  formulations 
with  special  emphasis  on  applications  to  metal  forming  processes.  The  literature 
dealing  with  specific  non-linearities  is  covered  in  Chapter  2  in  discussions  of  these 
non-linearities. 

The  rigid-plastic  finite  element  method  (FEM)  proposed  by  Lee  and  Kobayashi 
[12]  has  been  popularly  used  to  model  several  metal  forming  processes.  This 
method  uses  a  Lagrange  multiplier  to  enforce  the  incompressibility  condition  in 
pure  plastic  flow.  The  solution  scheme  in  the  method  is  based  on  using  an 
initial  guess  of  the  velocity  vector  with  small  perturbation  in  the  velocity  field 
in  subsequent  iterations  to  obtain  the  solution  as  a  function  of  time.  The  rate 
of  convergence  of  this  solution,  therefore,  depends  on  the  proximity  of  the  initial 
guess  to  the  exact  solution.  The  rigid-plastic  model  is  perhaps  locally  valid  for 
large  plastic  flow  where  the  solution  is  only  applied  to  regions  with  unrestricted 
plastic  flow.  Keeping  these  characteristics  of  the  rigid-plastic  method  in  mind, 
the  application  of  the  method  to  modelling  processes  such  as  cold  forging  [13,14], 
extrusion  [15],  strip-rolling  [9],  and  ring-compression  [  16]  is  quite  acceptable. 


The  rigid-viscoplastic  finite  element  method  originated  from  the  rigid-plastic 
method  due  to  the  work  of  Oh  wh">  later  used  it  to  model  forming  processes  with 
dies  of  general  shape  [17]  and  for  hot-forging  of  an  alloyed  engine  disk  (IS). 

Although  finite  strains  have  been  used  in  rigid-plastic  (or  viscoplastic)  analysis 
[19,20],  the  large  strain  formulations  have  been  more  commonly  used  with 
the  elastic-plastic  (or  viscoplastic)  material  model  by  a  number  of  researchers. 
Prominent  among  these  are  the  works  of  Win  [21]  and  Yamada  [22,23].  The  elasto- 
plastic  FEM  has  been  compared  to  the  rigid-plastic  FEM  by  Boer  [24]  for  upsetting 
of  cylinders.  Interestingly,  Boer’s  solutions  favor  the  rigid-plastic  FEM  for  the  4- 
node  isoparametric  elements.  However,  the  elastic-plastic  FEM  should  not  be 
condemned  on  this  evidence  alone,  since  the  upsetting  process  predominantly 
produces  areas  of  large  plastic  flow.  Secondly,  the  4-node  element  may  not  be 
the  proper  element  to  use  in  large  strain  applications,  where  the  “4CST-element” 
(a  rectangular  element  divided  into  four  constant  strain  triangular  elements  by 
joining  the  nodes  on  the  diagonals)  has  been  shown  to  have  a  more  accurate 
response  [25,26]. 

The  next  step  towards  a  more  accurate  modelling  of  metal  forming  processes 
involves  the  representation  of  friction  along  the  die-workpiece  interface.  Wifi  [21] 
accounted  for  friction  by  assuming  a  fully  sticking  contact  along  the  interface. 
Hartley  [27]  modelled  the  frictional  loading  condition  by  including  a  layer  of 
interfacial  elements.  However,  no  definite  criterion  for  establishing  the  properties 
of  these  elements  is  forwarded  in  the  paper.  More  realistic  models  of  the  frictional 
contact  conditions  are  used  by  Yamada  [22]  and  Kikuchi  [23].  The  latter  uses 
a  microscopic  definition  of  contact,  similar  to  the  one  proposed  by  Fredriksson 
[29].  However,  the  parameters  used  in  the  microscopic  model  do  not  seem  to  have 
significant  physical  properties  associated  with  them. 
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Thermal  stresses  have  also  been  coupled  with  large  plastic  deformations  by 
several  authors  to  analyze  forming  processes  such  as  upsetting  and  extrusion 
problems  [30,31].  However,  inspite  of  the  extensive  literature  on  analysis  of  forming 
processes,  no  paper  presents  a  convenient  large  strain  formulation  with  realistic 
representation  of  frictional  contact  in  one  package.  This  thesis  uses  a  combination 
of  conventional  large  deformation  elastic-plastic  formulation  with  a  simple  yet 
effective  frictional  contact  model. 

While  most  of  the  works  mentioned  previously  have  been  oriented  towards 
some  industrial  application,  several  authors  have  directed  their  attention  to 
developing  combined  non-linear  formulations  capable  of  both  analyzing  metal 
forming  processes  and  soiving  structural  mechanics  problems.  Two  popular 
approaches  have  emerged  from  the  efforts  of  Lee  [32]  and  Marcal  [33].  Lee 
proposed  the  theory  that  finite  strain  components  are  non-linear  expressions  in 
the  displacements  and  hence  may  not  be  additive  in  nature.  His  theory  using 
multiplicative  decomposition  of  the  deformation  gradient  has  been  followed  by 
Kleiber  [34,35!  among  others.  However,  the  additive  decomposition  of  the  strains 
into  elastic  and  plastic  parts  has  been  more  commonly  used  in  literature  and 
has  been  followed  in  this  thesis.  Following  Marcal’s  formulation  of  geometric 
and  material  non-linearities,  several  other  formulations  have  been  forwarded. 
Among  these  formulations,  the  popular  ones  are  due  to  Hibbit  [36],  Bathe  [37], 
McMeeking  [3S!,  and  Yamada  [39],  These  formulations  are  based  on  the  Prandtl- 
Reuss  equations  [40]  coupled  to  geometric  non-linearities  by  either  a  Lagrangian 
or  an  Euleria.n  formulation  using  energetically  conjugate  stress  and  strain  tensors. 
Relever.t  detaiis  of  these  formulations  are  given  with  references  in  Chapter  2. 

The  formulations  mentioned  in  the  preceding  paragraphs  can  be  used  to  model  a 
number  of  metal  forming  processes,  including  the  ausroiling  process.  However,  the 


unique  features  of  the  ausroliing  process  cannot  be  modelled  by  existing  software 
packages.  Hence,  using  the  existing  non-linear  formulations  as  a  starting  block, 
a  non-linear  finite  element  analysis  program  (NOFEAP)  has  been  developed  by 
the  author.  All  the  characteristics  of  ausroliing  have  not  been  modelled  in  the 
program.  This  thesis  describes  the  non-linear  formulations  and  the  model  of  the 
rolling  process  included  in  NOFEAP.  Various  benchmark  and  sample  problems 
solved  using  NOFEAP  are  also  included  in  the  thesis. 


Chapter  2 

NON-LINEAR  FORMULATIONS 


As  stated  in  the  previous  chapter,  the  ausrolling  process  has  several  characteris¬ 
tics  which  must  be  reflected  in  the  numerical  model.  Keeping  these  characteristics 
in  mind,  the  model  must  account  for  the  following:  (l)  non-linear  material  behavior 
due  to  plastic  flow;  (2)  non-linear  geometric  behavior  due  to  large  strains;  (3) 
frictional  contact  conditions  at  the  die-workpiece  interface;  (4)  roiling,  shearing 
and  swaging  operations;  (5)  temperature  effects  on  stress  distribution  and  material 
properties;  and  (6)  three-dimensional  analysis  capability.  Of  these  features,  the 
first  three  features  and  a  model  of  the  rolling  process  have  been  implemented  in 
the  computer  program  NOFEAP.  Shearing  and  swaging  operations,  which  are  part 
of  the  extension  to  three  dimensions,  and  the  thermal  effects  have  not -yet  been 
included  in  the  program. 

This  chapter  describes  the  non-linear  formulations  used  in  the  program. 
Description  of  each  form  of  non-linearity  is  preceded  by  a  review  of  the  relevent 
literature  on  the  subject.  The  finite  element  formulation  of  these  non-linearities, 
to  obtain  the  matrix  equations,  follows  their  description.  A  concise  narrative  of 
the  NOFEAP  computer  program  is  presented  at  the  end  of  this  chapter. 

2.1  Elastic-plastic  Material  Behavior 

In  the  following  sections,  discussion  of  formulations  and  relevent  references  is 
restricted  to  the  solution  of  plasticity  problems  by  the  finite  element  method.  The 
obvious  reason  for  such  a  limitation  is  that  the  earlier  techniques,  such  as  slip¬ 
line  theory,  mathematical  and  theoretical  approaches,  etc.,  have  been  completely 
overwhelmed  by  the  numerical  approach. 


2.1.1  Review  of  Literature 


Plasticity  is  one  of  the  fields  which  has  greatly  benefited  by  the  development  of 
the  finite  element  method.  Initial  approach  to  solving  the  plasticity  problems 
used  either  the  direct  stiffness,  or  the  initial  strain  method.  The  latter  was 
initially  formulated  using  the  finite  difference  method,  but  resulted  in  an  unstable 
solution.  It  was  later  adopted  by  Gallagher  [41]  to  the  finite  element  method. 
Simultaneously,  the.  tangent  stiffness  method  was  being  developed  by  Pope  [42]  and 
Marcal  [2]  to  solve  the  elastic-plastic  problem.  This  method  is  essentially  a  piece- 
wise  solution  to  the  non-linear  problem.  A  summary  of  these  early  developments 
in  numerical  solution  of  elastic-plastic  problems  has  been  given  by  Marcal  [43]  and 
Yamada  [44].  Yamada  [45]  also  introduced  a  plastic  stress-strain  matrix  obtained 
from  the  Prandtl-Reuss  relations.  The  use  of  this  and  other  matrix  formulations 
has  since  been  popular  in  the  non-linear  solutions.  Zienkiewicz  and  coworkers  [46- 
48]  introduced  the  “initial  stress*  approach  which  includes  a  general  development 
of  the  elastic-plastic  constitutive  tensor. 

The  formulation  used  in  this  thesis  is  similar  to  the  one  given  by  Zienkiewicz 
and  is  described  by  Owen  [5].  The  basic  assumptions  in  the  mathematical  theory 
of  plasticity  and  the  associated  expressions  are  outlined  in  the  following  sections. 
A  complete  treatment  of  the  mathematics  is  given  by  Hill  [49]. 

2.1.2  Mathematical  Formulation 

A  general  formulation  of  the  elastic-plastic  problem  is  given  in  this  section.  The 
problem  of  elastic-plastic  flow  is  characterized  by  a  yield  criterion  given  as 

/(<?)  =  K{k)  (la) 

W 


or, 


F{a,K)  =  f(a)  —  K{k)  =  0 
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The  right  hand  side  of  (la)  is  a  function  of  the  work  hardening  parameter  K. 
The  following  derivation  assumes  isotropic  expansion  of  the  yield  surface.  Similar 
derivation  for  kinematic  hardening  has  been  given  by  Dogui  [50j,  and  for  mixed 
hardening  by  Axelsson  [51,52]. 

The  elastic  behavior  of  the  material  is  represented  by  /  <  K ,  while  for  plastic 
flow  /  =  K.  At  a  plastic  state,  the  incremental  change  in  the  yield  function  due 
to  am  increment  in  stress  is  given  by  the  following: 


df  =  —da 
da 


(2) 


Then  if  df  <0,  elastic  unloading  occurs,  and  the  stress  point  returns  inside  the 
yield  surface;  if  df  =  0,  neutral  loading  occurs,  and  the  stress  point  remains  on 
the  current  yield  surface;  and  if  df  >  0,  plastic  loading  occurs  with  the  stress 
point  moving  to  the  envelope  of  an  expanding  yield  surface. 

Following  initial  yielding  of  material,  an  increment  of  strain  is  assumed  to  be 
the  sum  of  an  elastic  and  a  plastic  increment. 


de  =  de*  +  de} 


(3) 


A  contrasting  theory  for  large  strains,  based  on  the  decomposition  of  the 
deformation  gradient  into  elastic  and  plastic  parts,  has  been  forwarded  by  Lee 
[32,53].  A  discussion  of  this  controversy  is  also  provided  by  Lee  [54], 

In  order  to  derive  a  relationship  between  the  plastic  strains  and  the  stress 
increments,  the  plastic  strains  are  assumed  to  be  proportional  to  the  stress  gradient 
of  the  plastic  potential  Q. 


de p  =  A 


dQ 

da 


(4) 


Here  A  is  a  proportionality  constant  called  the  plastic  multiplier.  Using  the  yield 
function,  F,  for  the  plastic  potential  leads  to  the  associated  theory  of  plasticity 
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which  has  been  used  in  this  derivation.  In  such  a  case,  equation  (4)  can  be 
rewritten  as  follows: 

dF 

=  A-  (5) 

Defining  the  flow  vector,  a,  by 

_  dF 

a  =  jg  <6> 

a  differential  of  equation  (lb)  can  be  written  in  the  following  form: 


or, 

Here  A  is  defined  by 


dF  dF 

dF  =  —dc 7  +  ~z~dK  =  0 

dd  ok 


a?  dcr  —  AX  —  0 


A  =  - 

A  dK 


(7a) 

(7b) 

(8) 


For  the  von  Mises  criterion,  /  =  J2  the  second  deviatoric  stress  invariant,  and 
equation  (5)  reduces  to  the  familiar  Prandtl-Reuss  relations.  Since  the  elastic 
strain  increments  can  be  related  to  the  stress  increments  using  the  constitutive 
matrix  D ,  equation  (3)  can  be  rewritten  as  follows: 


dF 

de  =  [D)~ldd  +  A— 

OO 

Using  (9)  and  (7b),  the  following  expression  for  A  is  obtained: 

a  T[D\de 

A  — 


(9) 


(10) 


A-i-  a?\D)a 

Substutiting  (10)  in  (9)  yields  the  complete  elastic-plastic  stress-strain  relation  as 
given  next. 

dd  =  [Dep\de  (11) 

with 

\D<P\  =  [D\  -  (1=) 


A  +  ^[Dja 


For  the  work  hardening  hypothesis,  the  constant  A  is  seen  to  be  equal  to  the  slope 
[H  )  of  the  hardening  curve  [5],  The  elastic-plastic  constitutive  matrix,  given  by 
(12),  has  to  be  corrected  for  large  rotations  as  explained  in  the  following  sections. 

2.2  Geometrically  Non-linear  Behavior 

The  following  sections  review  the  formulations  defining  geometrically  non-linear 
behavior  and  present  the  total  and  the  updated  Lagrangian  formulations. 

2.2.1  Review  of  Literature 

Large  displacement  analysis  was  initially  of  interest  to  structural  engineers  in 
solving  bending  and  buckling  problems.  Marcai  [33]  presented  one  of  the  first 
formulations  combining  the  non-linear  material  behavior  with  large  displacements. 
He  has  also  given  a  historical  development  of  the  large  displacement  formulations. 
At  the  same  time  Lee  [32]  proposed  that  for  finite  elastic-plastic  strains,  the 
total  strain  may  not  be  the  sum  of  the  elastic  and  plastic  components  due  to 
the  dilatational  nature  of  the  strains.  However,  the  additive  decomposition  of 
the  finite  strains  is  still  popularly  used  and  has  been  followed  here.  Hibbit 
[36]  introduced  the  load  correction  terms  for  the  first  time  in  large  deformation 
analysis.  Hofmeister  [55]  presented  another  version  of  this  formulation  and  used  it 
to  analyze  large  displacements  of  structural  components.  Bathe  [37]  provided 
a  complete  description  of  the  total  Lagrangian  formulation,  which,  similar  to 
the  approach  of  Marcai  [33],  uses  the  initial  state  as  the  reference  configuration. 
Bathe  also  introduced  the  updated  Lagrangian  formulation  which  employs  the  last 
updated  position  as  the  reference.  Many  authors  erroneously  refer  to  the  updated 
Lagrangian  method  as  the  Eulerian  method.  The  subtle  difference  between  the  two 
is  explained  by  Gadala  [56].  Bathe  used  these  formulations  for  static  and  dynamic, 
large  deformation  analyses  of  structures  [57-59] .  Osias  [60]  and  McMeeking  [38] 
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forwarded  Eulerian  formulations  using  relative  velocities  as  independent  variables 
and  using  rate  equilibrium  equations. 


Tbe  rate  formulation  became  quite  popular  in  applications  involving  high  speed 
operations,  especially  metal  forming  processes.  Yamada  [39]  has  presented  a 
comprehensive  description  of  the  rate  formulation  with  explanation  of  various 
stress-rates  and  the  load  correction  matrices.  This  formulation  has  been  used 
to  analyze  large  flow  processes  such  as  extrusion  [23,39],  upsetting  and  stretch- 
forming  [21,22].  Several  authors  have  since  used  various  formulations  to  analyze 
different  metal  forming  processes  [28,37,61]. 

Although  a  number  of  formulations  have  been  forwarded  [36,28,39,59,62,63'., 
as  yet  there  is  no  universal  agreement  on  a  best  approach  to  solve  a  particular 
problem.  A  comprehensive  survey  of  the  available  formulations  is  given  by  Gadala 
[56,64].  Of  these  numerous  formulations,  the  total  and  updated  Lagrangian 
formulations  have  been  incorporated  in  the  program  NOFEAP.  Salient  features  of 
these  formulations  are  given  in  the  next  section,  the  details  being  available  in  the 
book  written  by  Bathe  [3j. 

2.2.2  Mathematical  Formulations 

In  the  analysis,  the  motion  of  a  body  is  considered  in  a  stationary  Cartesian 
coordinate  system,  as  shown  in  Fig  2.1.  All  kinematic  and  static  variables  are 
measured  in  this  coordinate  system.  Tensor  notation  is  followed  in  this  section. 

The  basic  equation  to  be  solved  is  given  by  (13),  which  expresses  the  equilibrium 
of  a  general  body  at  time  t  +  dt. 

f2a„  52E,,  W  =  -X  (13) 
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Here  <T ,y  and  EXJ  denote  the  engineering  stress  and  strain  tensors  respectively. 
With  the  standard  tensor  notation,  a  left  subscript  is  used  to  indicate  the  reference 
coordinate  axes  and  a  left  superscript  to  denote  the  configuration  of  the  body.  If 
a  quantity  occurs  in  the  same  configuration  in  which  it  is  also  measured,  the  left 
subscript  is  dropped.  Also,  an  incremental  quantity  is  denoted  by  dropping  the 
left  superscript.  For  instance,  the  term  q ux  indicates  the  itfl  component  of  the 
displacement  vector  at  time  t  4-  dt,  measured  with  respect  to  the  configuration 
at  time  0;  the  quantity  1Ut-  denotes  the  component  at  time  t  with  respect  to  an 
updated  configuration,  also  at  time  f;  and  oU^  indicates  incremental  displacement 
component  at  the  current  state  measured  with  respect  to  the  initial  state. 

The  following  sections  describe  the  derivation  of  the  final  form  of  equilibrium 
equations  for  the  total  Lagrangian  formulation.  The  equations  for  the  updated 
Lagrangian  formulation  are  presented  without  proof  since  these  can  be  obtained 
simply  by  following  the  steps  given  in  the  next  section. 


2.2.2. 1  Total  Lagrangian  Formulation 


The  stress  and  strain  tensors  selected  for  the  total  Lagrangian  (TL)  formulation 
are  the  2nd  Piola-Kirchoff  stress  tensor,  Sxj,  and  the  Green-Lagrange  strain 
tensor,  6ty,  respectively.  The  equilibrium  equation  (13)  takes  the  following  form 
using  these  tensors  with  the  initial  configuration  used  as  the  reference. 

f  IS ,,  Sic,,  °dV  =  (H) 

°v 


In  equation  (14),  the  Green-Lagrange  strains  are  uefined  by  the  following: 


G,  =  - 


o'-*; 


o“«.> 


ou;,» 


ouk,i  *  o uk,j 


(15) 


where, 


l 

o 


u 


(16) 


a 


■.V 
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The  stresses  and  strains 

are  decomposed  as  follows: 

o'S'.y  =  o  ^  ij  +  o  Sij 

(17a) 

and 

0£»y  =  offj  ~h  o^y 

(176) 

The  incremental  strains  are  further  assumed  to  be  the  sum  of  linear  and  non-linear 
parts  as  given  next: 


0 £\j  ~  0eij  "h  oVij 

(18) 

where, 

* 

0«i/  ~  ~(oui,3  d*  0 uj,i  +  ouk,x  *  o^k^j  +  o ^k,i  *  ouk,i ^ 

(19a) 

and 

1 

(196) 

oVij  ~  *  0 uk,j 

Equation  (14)  can  now  be  expanded  to  the  following  form: 


1 0Si,  <50e<y  °dV  +  J  JS,-,-  60rhi  °<iV  = 

o  y  o  y 

~  J  l0Sij  60eij  °dV  (20) 

0  V 

Equation  (20)  is  used  with  a  constitutive  law  and  is  converted  to  a  matrix  form 
to  be  solved  iteratively.  The  external  virtual  work,  2  R,  is  given  by  the  following: 

7  =  /  7*  Sui  °dV+  f  7*  6ui  °dA  (21) 

o  y  o 

Here  and  /j*  denote  the  components  of  body  forces  and  surface  tractions 
respectively. 

2. 2. 2. 2  Updated  Lagrangian  Formulation 

In  this  formulation,  the  Cauchy  stress  tensor,  TXJ,  is  used  with  the  following 
strain  definition: 

l^i j  ~  =  leij  "h  itfij  [2~d) 


>.V.V 


\  i  ,1  *  K 


t  vVv.VvP’/.'.'.V.'.'V.v /  v 
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where, 


i e«7  “  2^lU,,J  +  lUjr*’^ 


(226) 


The  equilibrium  equation  (13)  is  written  as  follows,  using  the  configuration  ai  time 


t  as  reference. 


I  ISij  8\c,j  'dV  =  2JJ 


Using  the  steps  outlined  in  the  previous  section,  equation  (23)  can  be  converted 
to  the  following  form: 

J  iSij  Slei3-  1dV  +  J  lT{j  Stfij  xdV  = 

11/  IV 

I'-,;  5 :e,'y  W  (24) 


The  external  virtual  work  is  given  by  the  following: 


'*  =  J  ■}■  Sui'dV+  hj!  SUi\ 


In  equation  (20),  05,y  is  replaced  by  {0 Dijrs  o^rs}  a^d  *m  (24),  iSij  is 
replaced  by  {x.D,yr3  lfira}  using  appropriate  incremental  constitutive  tensor 
Dxjr3  at  time  t  referred  to  the  initial  and  updated  configurations  respectively. 
The  next  section  describes  the  relation  of  the  constitutive  tensor  to  the  elastic  (or 
elastic-plastic)  constitutive  matrix  described  in  section  2.1. 

2.2.3  Use  of  Constitutive  Relations 


In  order  to  obtain  a  sensible  solution  to  a  given  problem,  it  is  imperative  that 
the  employed  formulation  contains  appropriate  model  of  both  the  kinematic  and 
constitutive  description.  The  kinematics  of  the  formulations  have  been  discussed 
in  the  previous  sections.  This  section  outlines  the  use  of  constitutive  relations  [3] 
for  the  total  and  the  updated  Lagrangian  formulations. 
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2.2.3. 1  Elastoplasticity  With  Total  Lagrangian  Formulation 

The  kinematic  description  of  this  formulation,  as  described  in  section  2. 2. 2.1, 
permits  large  displacements,  rotations  and  strains.  However,  the  formulation  is 
most  effective  for  large  displacements  and  rotations,  but  small  strain  analysis  of 
elastic-plastic  materials  [3].  Such  a  model  is  directly  obtained  by  substituting  the 
2  Piola-Kirchoff  stresses  and  Green-Lagrange  strains  for  engineering  stresses 
and  strains  in  the  elastic-plastic  formulation.  The  yield  function  is  written  as 


^(X-,  1K)=0 


(26) 


and  the  flow  rule  relates  an  increment  in  the  Green-Lagrange  strains  to  the  2 
Piola-Kirchoff  stress-gradient  of  the  yield  function. 

d'F 


nd 


dev-  =  l\ 

d'S. 


(27) 


The  stress-strain  relation  takes  the  following  form: 


dSij  =  0Dijra(der3  -  de.vT3) 


(28) 


and  an  expression  for  elastic-plastic  constitutive  matrix,  similar  to  equation  (12), 
is  obtained. 

2. 2. 3. 2  Elastoplasticity  With  Updated  Lagrangian  Formulation 

For  large  displacement,  rotation  and  strain  analysis  of  elastic-plastic  materials, 
the  updated  Lagrangian  method  with  Jaumann  stress  rate,  f,y ,  (ULJ)  formulation 
has  been  used  in  the  program  NOFEAP.  The  yield  function  is  now  a  function  of 
Cauchy  stresses  and  can  be  expressed  as  follows: 


1F(1rij,IK)  =  0 


(29) 
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The  flow  rule  relates  the  plastic  strain  increment  to  the  Cauchy  stress  derivative 


of  the  yield  function. 


dlF 


The  stress-strain  equation  (31)  relates  an  increment  in  the  Jaumann  stresses  to 
an  elastic  strain  increment. 

fijdt  —  iDijra{dsr3  de?3)  (31) 


In  (31),  d&ij  are  increments  of  the  deformation  tensor  and  are  the  linear  compo¬ 
nents  of  the  Euler  or  Almansi  strain  tensor  [65 j.  The  Cauchy  stress  increments 
are  computed  from  the  Jaumann  stresses  using  the  following  expression: 

X  —  Tijdt  ~r  /tp  Clpjdt  t  ~jp  ftpidt  (3 2) 

where  1Clij  a~re  components  of  the  spin  tensor. 

,  1  dlilj  dlU: 

n<i  =  2 (3uT  _  (33) 

Equations  (31-33)  are  used  to  compute  the  Cauchy  stresses  at  each  increment. 
In  order  to  obtain  the  constitutive  matrix  relating  the  2nd  Piola-Kirchoff  stress 
increment,  iS{j,  to  an  increment  of  deformation  tensor,  dC{j,  a  relation  between 
the  2nd  Piola-Kirchoff  stresses  and  the  Truesdell  stresses,  dalnn,  is  used  [36 j . 

,,  d°Xi  d°Xj  , 

dSii  =  wd^Zd^;da'mn  (34) 

where  liX]  is  the  determinant  of  the  deformation  gradient.  The  Truesdell  stress 


increments  are  related  to  the  Jaumann  stresses  by  the  following  relation: 

Tijdt  —  dcrmTl  -r  TrnpdCpn  t  Tnpd6pm  Trnnd&pp 


Equation  (31),  (34)  and  (35)  provide  the  desired  constitutive  relation.  The  term 
~mndepp  contributes  an  unsymmetric  matrix  to  the  stiffness  matrix.  However, 


depp  denotes  the  volumetric  strain  which  is  small  for  large  plastic  flow.  Also, 
the  ratio  of  lTmn  to  the  elastic  modulus  is  within  the  order  of  elastic  strain 
and  is  generally  negligible  in  finite  strain  analysis.  Therefore,  this  term  may  be 
dropped  from  the  relation  without  significantly  affecting  the  solution  while  saving 
considerable  computational  effort. 

The  external  virtual  work  contains  contributions  due  to  body  and  surface  forces, 
as  given  in  equations  (2 1)  and  (25).  Some  of  the  surface  loads  may  be  applied  by 
a  die  penetrating  the  workpiece.  A  treatment  of  such  contact  loads  is  given  in  the 
next  section. 

2.3  Surface  Non-linearities 

The  next  section  includes  a  review  of  the  formulations  used  to  numerically  solve 
a  contact  problem.  A  discussion  of  the  formulation  used  in  NOFEAP  follows. 

2.3.1  Review  of  Literature 

Although  the  analysis  of  contact  problems  has  been  an  area  of  investigation 
since  Signorini  [66]  first  applied  a  variational  method  to  it  in  1959,  numerical 
methods  were  not  actively  used  in  the  field  for  more  than  a  decade.  The 
finite  element  method  was  introduced  into  the  realm  of  contact  problems  by 
Chan  [4 1  in  1971.  Chan  interpolated  the  stiffness  at  the  point  of  contact, 
which  was  permitted  to  be  between  any  two  nodes,  from  the  stiffnesses  at  these 
nodes.  The  procedure  was  developed  for  linear  triangular  elements.  Similar 
approach  for  higher  order  elements  seems  complicated  and  has  not  been  followed 
in  subsequent  research.  Since  1971,  a  number  of  formulations  and  algorithms 
[6,29,67-71!  have  been  presented.  Most  of  these  formulations  use  interface  or 
bond  elements  on  the  contact  surface.  An  entirely  different  approach  has  been 
taken  by  Fredriksson  [29,72!,  who  developed  a  formulation  similar  to  the  theory 
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of  plasticity.  He  introduced  a  yield  criterion  and  a  flow  rule  into  contact  problem 
definition.  This  formulation  has  been  adopted  by  Kikuchi  [28)  and  Cheng  (31). 
Both  Fredriksson  [73]  and  Cheng  [3l]  have  presented  experimental  means  of 
determining  the  frictional  properties.  However,  some  of  the  parameters  used  in 
this  micromechanical  formulation  seem  to  lack  physical  significance. 


These  algorithms  have  been  used  to  analyze  a  variety  of  problems,  ranging 
from  indentation  [29]  and  contact  in  machine  elements  [72,74]  to  dynamic  contact 
[70,75,76].  The  method  used  in  this  thesis  is  developed  by  Okamoto  [6]  and  is 
described  next. 


2.3.2  Mathematical  Formulation 


In  this  formulation  [6],  a  contact  element  is  used  at  the  surface.  The  stiffness 
matrix  and  load  vector  of  the  element  are  generated  using  the  equilibrium 
equations  and  continuity  conditions  imposed  on  the  force  and  displacement 
increments,  as  given  next. 

Let  uf’,  Rf’  and  uf ,  R ^  denote  the  displacements  and  reactions  at  the 
contact  surfaces  of  bodies  A  and  B  respectively,  defined  along  the  local  normal 
and  tangential  coordinates  (Fig  2.2).  A  prefixed  A  denotes  an  increment  in  the 
quantity.  The  constraints  for  various  contact  conditions  are  the  following: 

1.  Open  contact  condition: 


A  R?  =  -A  R?  =  0 
A  uf  -  A  uf  =  A  U 

2.  Sticking  contact  condition: 

A  R?  =  -  A  R?  =  A  Ri 
A u*  -  A =  <5n 
A  uf  —  Auf  =  0 


(36a) 

(36b) 

(37a) 

(376) 

(37C) 
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3.  Sliding  contact  condition: 

A  R£  =  -&Rn  =  A  Rn 

(38a) 

A  R?  =  -A  R?  =  ±/zA  Rn 

(386) 

A u*  -  Au%  =  Sn 

(38c) 

Auf  —  Attf  =  Alt 

(38<f) 

The  nodal  unknowns  are  the  incremental  displacements  and  the  relative 
increments  A whenever  the  motion  along  the  direction  is  permitted,  or  the 
incremental  contact  forces,  if  no  such  motion  is  allowed.  Thus,  although  there  are 
six  variables  associated  with  each  contact  element,  compared  to  four  in  commonly 
used  elements,  the  nodal  reactions  and  relative  increments  are  obtained  as  part  of 
the  solution. 

The  nodal  displacement  increments  All,-  and  the  contact  force  increments  A i?t 
under  the  known  contact  conditions  Cn  and  a  test  load  increment  A F  are  obtained 
by  solving  the  combined  matrix  equation.  This  system  of  simultaneous  equations 
is  obtained  by  assembling  the  structurak^tiffness  matrices  and  load  vectors  with 
the  contact  stiffness  matrices  and  load  vectors.  The  solution  vectors  are  used  to 
calculate  the  load  increment  ratio  a k  at  the  kth  step  as  given  by 


and 


A  Fi  =  afc  *  A  .Ft  (39a) 

Au,-  =  Qfc  *  An,-  (396) 


The  factor  Qfc  is  the  smallest  fraction  of  the  load  required  to  cause  a  change 
in  the  contact  statuses  of  one  contact  element.  A  discussion  of  the  load  factor 
computation  scheme  is  given  by  Okamoto  (6l,  whose  work  has  been  extended  by 
Saiamon  (77).  Expressions  for  these  load  factors  are  presented  in  Appendix  A. 
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2.4  Finite  Element  Formulation 


The  equations  (20)  and  (24)  are  coupled  with  elastic-plastic  constitutive  laws 
and  constraint  equations  at  contact  surface.  Both  in  elastic-plastic  solution 
and  large  deformation  analysis,  the  out  of  balance  virtual  work  is  reduced  by 
performing  iterations  until  the  difference  between  the  external  and  internal  virtual 
work  is  negligible  within  a  certain  convergence  criterion.  For  the  total  Lagrangian 
formulation,  equation  (20)  is  rewritten  as  follows  and  is  solved  for  an  iteration 
step  k. 


1 0 Dijr,  „e<?  <0 e.-,-  °dV  +  j  ISij  S0T}ij  °dV  = 


°V 


-  !  is!} 


°v 

(k-l)  c  Ak- 1)  0 


dV 


(40) 


Equation  (24)  for  the  updated  Lagrangian  formulation  is  written  in- the  following 
form: 


h 


D 


IJT3 


Jk) 
1  ers 

8xeij  1dV  -f- 

I'm  « 

> V 

■z  - 

/ 

7y[k-i) 

ldV 


(41) 


where  the  constitutive  matrices  o-^ij'rs  and  i Dijrs  are  obtained  using  the 


method  outlined  in  section  2.2.3.  The  following  development  is  only  for  the  total 
Lagrangian  formulation.  Similar  expressions  can  be  obtained  for  the  updated 
Lagrangian  method  by  following  these  steps.  Using  the  standard  finite  element 
discretization  technique  [7 j ,  equation  (40)  can  be  converted  to  the  following  matrix 
equation: 

\[Kl\  -  \KsA\  (i«)  =  {o-R}  -  {oR'*"’}  («“) 


»■»  ■ 


m 


m 


p 


m 


■»  A 


m 

m 


[[■K’l]  +  [Knl\]  \Au}  =  {2Ai?(A:)} 


(426) 


where, 

[ Kl]  =  linear  stiffness  matrix  —  [i^s]  + 

\Ks\  —  small  displacement  stiffness  matrix; 

[ifs]{Au}<5{Au}  =  J  0e,a  0 Dijrs  S0e^  °dV  .  (43) 

°V 

[Kq\  =  geometric  stiffness  matrix; 


[KcI{Au}S{au}  =  J  04,  „Dijr,  S0e%L  °dV  + 

o  V 

J  0e?,L  ,Dijr,  60e%  °dV  +  j  „c-vsL  0Dijr,  Sae"L  °dV  (44) 

o  y  o  y 

[K nl\  =  non-linear  (initial  stress)  stiffness  matrix; 


[^^l]{Au}6{Au}  =  J  lSij8QTjij°dV  (45) 

°V 


{qJ2}  =  total  external  virtual  work  at  time  t  ~  dt; 

=  total  internal  work  at  time  t  +  dt,  but  at  previous  iteration  step; 
=  increment  in  external  virtual  work  at  time  t  +  dt  and  at  kth  step. 
In  (43)  and  (44),  the  quantities  and  are  defined  as  follows: 


osf'j  -  ~(oui,j  +  o uj,i)  (46a) 

0ei ij  =  “(o ufc,»  *  0 uk,j  T  o uk,i  *  oUk>j)  (466) 

z 

The  non-linear  strain  increments,  o7i;  >  are  defined  by  equation  (19b).  The 
expressions  for  updated  Lagrangian  formulation  are  similar  to  the  ones  given  here 
except  for  the  absence  of  the  geometric  stiffness  matrix.  Equation  (42b)  can  be 
expressed  as 


(A']{au}  =  =  {AiJ} 


(47) 
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with  the  Incremental  virtual  work  {aA}  given  by  the  following: 


{aA}«5{au}  =  f  A fb8ui°dV  +  J  Af*5ui°dA 


If  the  contact  surface  of  the  body  is  separated  from  the  remaining  surface,  the 


incremental  virtual  work  may  be  broken  up  as  follows: 


{ar}  =  {Air6}  +  {Air} 


Here  {Ait!c}  contains  the  work  of  the  unknown  contact  forces,  and  {aAs6} 


includes  the  contributions  of  the  body  forces  and  known  surface  tractions.  For 


two  contacting  bodies,  A  and  B,  equations  (47)  and  (49)  are  written  as 


=  {A -r  {aJ%} 


(50a) 


[Kb\{  A*b}  =  {A  R%b}  +  {AiJJ,} 


(506) 


{Ail|}  =  -{A^} 


(50c) 


Equations  (50)  can  be  combined  to  yield  a  single  matrix  equation  as  follows  [6]: 


0  Kb  KqCc  J  I  Aub  I  =  I  a R3Bb 
Kca  A cb  A ccJ\aQ  j  \  6c 


Aua\  fAR3Ab 


where,  KACi  A BC,  KcAi  Kcb>  and  A cc  are  contact  stiffness  matrices 

_  / 

obtained  from  constraint  equations  (36-38).  A Q  are  incremental  vectors  of 
contact  forces  or  relative  displacements  and  Sc  are  vectors  of  initial  gaps  between 


nodal  pairs  along  the  contact  surface. 


Equation  (51)  includes  non-linearities  due  to  large  deformations  in  the  structural 
stiffness  matrices,  KA  and  K b\  due  to  contact  conditions  in  establishing  the 


contact  statuses  at  each  step;  and  due  to  plastic  flow  in  the  stiffness  and  stress 


I 


u«  v  ■  wwwmi  wwiuifw  mngmgrnsiia] 


22 

computations.  This  can  be  expected  to  lead  to  convergence  problems.  Programs 
with  conventional  algorithms  may  lack  the  flexibility  required  to  solve  such 
problems.  A  macro  programming  language  [7]  overcomes  this  difficulty  by  letting 
the  user  control  the  solution  algorithm.  A  brief  description  of  this  feature  is  given 
in  the  next  section. 

2.5  Program  Features 

The  non-linear  formulations  presented  in  previous  sections  have  been  used  to 
develop  a  non-linear  finite  element  analysis  program  (NOFEAP).  The  program 
has  been  designed  as  a  flexible  tool  for  research,  educational,  or  applications 
environments,  which  demand  frequent  modifications  to  meet  each  new  problem’s 
requirements.  The  program  has  been  developed  on  a  VAX  11/7S2  running  under 
VMS  version  4.3  operating  system,  and  is  also  operational  on  a  Data  General 
Eclipse  MV/10000  with  AOS/VS  7.54  operating  system.  The  program  consists  of 
several  general  modules:  (l)  Problem  control;  (2)  Problem  definition  and  mesh 
input;  (3)  Element  library;  (4)  Problem  solution;  and  (5)  Graphics  outputs. 

The  problem  solution  module  is  centered  around  a  unique  macro  programming 
language  [7]  concept  in  which  the  solution  algorithm  is  written  by  the  user.  There 
are  sufficient  macro  instructions  included  in  the  system  for  many  applications  in 
structural  mechanics;  however,  each  user  may  modify  the  program  by  adding  new 
language  features  to  meet  the  user’s  specific  application’s  requirements. 

The  NOFEAP  system  includes  a  general  element  library.  Elements  are  currently 
available  to  model  linear  and  non-linear  problems  in  structural  mechanics. 
Included  in  the  existing  library  are  elements  for  solving  linear  elasticity  problems 
in  plane  stress,  plane  strain,  or  axisymmetric  analysis.  Non-linear  elements  are 
available  for  an  elastic-plastic  solution  of  planar  problems  using  a  variety  of  yield 
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criteria.  The  total  and  the  updated  Lagrangian  formulations  are  modelled  in 
elements  for  geometrically  non-linear  analysis  of  planar  elements.  Sun  ace  non- 
linearities  have  been  modelled  by  including  the  contact  problem  formulation 
described  in  the  previous  section. 

External  loads  can  be  applied  in  the  problems  by  prescribing  the  forces  or 
displacements  as  well  as  by  rigid  body  translation  or  rotation  of  contacting  bodies. 
The  rolling  process  has  not  been  modelled  in  any  existing  commercial  program, 
although  some  programs  do  approximate  the  rolling  process  by  assuming  the 
pressure  distributions  along  the  contact  surface  [10,1 1].  The  features  of  NOFEAP 
and  descriptions  of  the  available  macro  commands  are  given  in  the  NOFEAP 
Users’  Manual  [8j. 

The  program  NOFEAP  has  been  used  to  analyze  the  deformations  in  a  cylinder 
when  it  is  rolled  against  a  hard  cylinder.  The  plastic  flow  and  stress  distribution 
in  a  gear  tooth,  when  the  gear  is  roiled  against  a  rigid  die,  has  also  been  studied. 
These  results  are  presented  in  Chapter  4.  The  implementation  of  each  of  the 
non-linear  formulations  in  NOFEAP  has  been  verified  by  solving  a  number  of 
benchmark  problems  -  problems  which  have  been  solved  both  analytically  and 
numerically  by  various  authors.  Some  of  the  problems  have  also  been  solved 
experimentally  by  different  authors.  The  results  of  some  of  these  verification 
exercises  are  presented  in  the  next  chapter. 
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Chapter  3 

ANALYSIS  OF  SOME  NON-LINEAR  PROBLEMS 


The  non-linear  formulations  in  NOFEAP  have  been  checked  for  accuracy  by 
solving  a  number  of  benchmark  examples.  These  problems  have  been  selected  to 
test  the  formulations  both  individually  and  in  various  combinations.  Results  of 
some  of  these  exercises  are  given  in  this  chapter. 

The  program  NOFEAP  is  capable  of  using  a  number  of  planar  isoparametric 
elements,  namely,  three  through  nine  node  elements.  The  analysis  can  be  carried 
out  for  plane  strain,  plane  stress  or  axisymnetric  bodies.  The  elastic-plastic 
analysis  can  be  performed  for  materials  obeying  either  Tresca,  von  Mises,  Drucker- 
Prager  or  Mohr-Coulomb  criterion  for  yielding.  However,  the  results  are  only 
included  for  von  Mises  material  since  this  model  is  to  be  used  to  study  the  gear- 
rolling  process. 

3.1  Thick  Cylinder  Under  Internal  Pressure 

The  first  example,  analysis  of  a  thick  cylinder  subjected  to  an  internal  pressure, 
provides  a  simple  means  of  comparing  the  numerical  results  with  analytical 
ones.  This  problem  has  been  solved  by  several  authors  [5,47,78]  and  the  stress 
distributions  are  also  available  in  the  literature. 

The  finite  element  mesh  used  in  the  analysis,  consisting  of  13  eight-node 
elements,  is  shown  in  Figure  3.1.  Only  a  quadrant  of  the  cross-section  needs  to  be 
modelled  due  to  symmetry  about  both  coordinate  axes.  The  material  properties 
are  assumed  to  be  the  following: 


modulus  of  elasticity,  E 
yield  s  iii  CSS,  y 


=  0.210E  -  06 N/mmr\ 

=  240.0 OjV/mm:; 


UNDER  INTERNAL  PRESSURE 


slope  of  hardening  curve,  H  (=A)  =  0.00 N/mm2; 

Poisson’s  ratio,  v  =  0.30. 

Figure  3.1  also  displays  the  number  of  nodes  and  elements  along  with  the  number 
of  integration  points  (G.P.)  used  in  the  problem.  The  cylinder  is  subjected  to  an 
internal  pressure  of  80 N/mm7,  which  is  then  doubled  in  steps  of  20N  /  mm7 . 
Non-dimensionalized  curves  obtained  from  the  analysis  are  presented. in  Fig  3.2- 
3.4.  The  computed  stresses  referred  to  in  this  thesis  imply  the  stresses  at  Gauss 
points  unless  otherwise  noted.  The  distances  and  displacements  are  normalized 
by  the  inner  radius,  r(a),  while  the  pressure  and  stresses  are  normalized  by  the 
yield  stress,  <Xy.  The  shear  modulus  of  the  material  is  denoted  by  G. 

The  loading  curve  of  the  problem  is  given  in  Fig  3.2  along  with  the  numerical 
solution  of  Nayak  [47]  and  the  analytical  solution  of  Hodge  [78].  The  solution 
obtained  by  NOFEAP  exactly  matches  the  analytical  solution  while  differing 
slightly  from  the  numerical  results  of  Nayak.  The  slight  variation  between  the  two 
numerical  solutions  may  be  attributed  to  a  difference  in  convergence  tolerance  or 
criterion,  numerical  formulation  and  finite  element  mesh.  However,  the  agreement 
of  NOFEAP  solution  with  the  analytical  results  indicate  the  accuracy  of  the 
elastic-plastic  formulation. 

The  hoop  stress  distribution  across  the  thickness  of  the  cylinder  is  presented 
in  Fig  3.3.  The  results  of  NOFEAP  compare  exceedingly  weil  with  the  analytical 
ones  while  again  deviating  from  Nayak’s  solutions.  Finally,  the  stress  distributions 
for  various  steps  are  presented  in  Fig  3.4.  The  behavior  of  the  curves  is  as  expected 
with  the  maximum  hoop  stress  shifting  away  from  the  inner  surface  with  increasing 
plastic  flow  while  the  radial  stresses  maintain  a  traction  free  outer  surface. 
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PRESSURE  VERSUS  DISPLACEMENT  CURVE  FOR  THICK  CYLINDER 


3.2  Small  Strain  Extrusion 


This  example  demonstrates  the  application  of  the  elastic-plastic  analysis  to  a 
practical  problem:  the  extrusion  process.  The  actual  process  has  been  simplified 
by  assuming  a  frictionless  die  and  small  plastic  strains  so  that  large  deformation 
formulations  and  frictional  contact  conditions  do  not  have  to  be  included  in  the 
analysis.  This  problem  has  also  been  solved  by  Nayak  and  Zienkiewicz  [47]  who 
have  compared  their  solutions  to  an  experimental  study. 

The  finite  element  mesh  (Fig  3.5)  employs  four-node  plane  strain  elements  in  the 
analysis.  Symmetry  of  the  die  about  the  vertical  axis  in  the  figure  permits  the  use 
of  half  the  section.  An  elastic,  but  perfectly  plastic  material  is  considered  in  the 
extrusion  process  with  an  elastic  modulus  of  68,950  MPa,  a  yield  stress  of  68.95 
MPa  and  a  Poissons  ratio  of  0.3.  Forward  extrusion  is  carried  out  by  steadily 
forcing  the  die  into  the  wider  end.  Since  only  small  strain  analysis  is  performed; 
the  final  displacement  of  the  die,  approximately  0.762  mm,  is  negligible  compared 
to  that  in  a  typical  process.  However,  the  analysis  does  produce  results  which  can 
be  compared  to  both  numerical  and  experimental  results. 

The  contours  of  effective  stress  are  shown  in  Fig  3.6.  These  contours  are 
indicative  of  the  progress  in  the  deformation  of  the  material  since  they  point 
to  the  origin  and  the  direction  of  plastic  flow.  These  contours  do  agree  with  those 
of  Nayak  [47!  which  could  not  be  reproduced  here.  The  largest  amount  of  flow 
(Fig  3.7)  occurs  at  the  corners  near  the  bend  in  the  die,  as  can  be  expected.  The 
load  versus  deflection  curve  of  NOFEAP  analysis  has  been  compared  to  that  of 
Nayak  in  Fig  3.8.  The  pressure,  P,  is  computed  by  summing  the  reactions  at  the 
penetrating  punch  and  dividing  the  sum  by  the  width  at  the  die.  As  can  be  seen 
in  Fig  3.8,  the  two  results  are  in  excelled  agreement. 
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Vcsh  for  Finite  Element  Analysis 
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ELASTIC  -  PLASTIC  FORWARD  EXTRUSION 


3.3  Large  Displacement  Analysis  of  a  Cantilever 

There  is  a  dual  purpose  in  selecting  large  displacement  analysis  of  a  cantilever 
to  test  the  accuracy  of  the  geometrically  non-linear  formulations.  The  problem 
has  been  solved  both  analytically  (Holden  [79})  and  numerically  (Bathe  [59]). 
It  not  only  provides  a  comparison  between  these  results  and  those  obtained 
by  the  total  Lagrangian  (TL)  and  the  updated  Lagrangian  (ULJ)  formulations, 
but  also  between  the  TL  and  ULJ  methods  themselves.  The  latter  comparison 
is  significant  since  with  proper  implementation  of  the  two  formulations  with 
appropriate  constitutive  relations,  identical  numerical  results  should  be  obtained 
using  either  formulations. 

The  cantilever  is  modelled  using  eight-node  isoparametric  elements.  The  mesh 
and  material  properties  are  given  in  Fig  3.9.  The  beam  is  254  mm  long,  and  25.4 
mm  In  width  and  depth.  The  elastic  modulus  of  the  material  is  82.74  MPa  and 
the  Poisson’s  ratio  is  0.2.  The  yield  stress  is  set  to  a  sufficiently  high  value  to 
maintain  elastic  response  in  the  material. 

The  beam  is  subjected  to  uniformly  distributed  loads  along  the  top  and  the 
bottom  surfaces.  The  pressure  is  increased  in  steps  of  689.5  Pa  for  100  increments. 
The  displacement  of  the  beam  at  13,790  Pa  is  shown  in  Fig  3.10.  A  comparison 
of  the  various  loading  curves  is  given  in  Fig  3.11,  which  shows  the  variation  in  the 
vertical  displacement  at  the  tip  with  increasing  pressure. 

A  linear  (small  displacement)  solution  is  superimposed  on  other  curves  in 
Fig  3.11  to  emphasize  the  difference  in  geometrically  non-linear  analysis  from 
a  linear  analysis.  It  should  be  noted  that,  the  results  of  NOFEAP  for  TL  and 
ULJ  formulations,  with  the  load  applied  in  100  steps,  match  exactly  with  each 
other,  indicating  a  consistent  implementation  of  these  formulations  in  the  program. 
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These  results  compare  extremely  well  with  those  of  Bathe  [59],  The  irregularity  in 
the  reproduced  curve  of  Bathe’s  is  ascribed  to  human  error  in  reading  the  curve 
from  the  published  article.  The  analytical  results  of  Holden  [79]  are  not  reproduced 
here,  since  these  have  been  shown  [59]  to  match  Bathe’s  solutions. 

The  agreement  of  the  non-linear  formulations  of  NOFEAP  and  the  existing 
analytical  as  well  as  numerical  solutions  is  found  to  be  satisfactory.  As  a  further 
test,  the  load  was  applied  in  only  five  steps  and  convergence  iterations  were 
performed  for  TL  formulation.  The  load-deflection  curve  obtained  for  this  case 
(Fig  3.11)  still  compares  fairly  well  with  other  results,  indicating  further  the 
accuracy  and  reliability  of  the  formulation.  In  the  next  example,  both  the  material 
and  geometric  non-linearities  are  combined  in  one  problem. 

3.4  The  Upsetting  Process 

The  upsetting  of  a  circular  cylinder  is  characterized  by  large  plastic  deforma¬ 
tions.  As  a  benchmark  problem  for  testing  the  validity  of  numerical  methods 
for  analysis  of  metal  forming  process,  it  has  been  solved  by  researchers  both 
numerically  [22, 31]  and  experimentally  [24] . 

The  cylinder  is  initially  20  mm  in  radius  and  30  mm  in  height.  The  mechanical 
properties  of  the  material  are  as  follows:  elastic  modulus  =  200  GPa,  Poisson’s 
ratio  =  0.3,  yield  stress  =  700  MPa,  and  slope  of  hardening  curve  —  300  MPa. 
Only  a  quarter  of  the  cylinder  is  analyzed  using  “4-CST”  elements  (Fig  3.12)  due 
to  symmetry.  Fully  sticking  conditions  are  assumed  from  the  top  of  the  block. 
The  block  is  deformed  by  compressing  at  the  top  in  steps  of  0.25%  of  initial  height 
to  reach  a  final  reduction  of  30%  in  height.  A  second  study  using  0.50%  reduction 
in  height  as  step  size  has  also  been  conducted  to  observe  the  effect  of  load  step  on 
the  stresses  and  loads. 
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The  deformed  cylinder,  Fig  3.13,  shows  the  bulging  at  the  centerline,  as  has  been 
reported  in  the  experiments  of  Beer  [2-4] .  The  loading  curves  for  the  upsetting 


process  for  various  cases  are  shown  in  Fig  3.14,  along  with  curves  for  two  different 
load  steps  computed  by  Yamada  [22’ .  The  maximum  error  between  Yamada’s 
solution  and  NOFEAP  solution  for  a  step  size  of  0.25 %  has  been  found  to  be 
1.2~  -  an  acceptable  value  considering  the  non-linearities  involved  in  the  problem. 
Although  Yamada  [22  (and  Cheng  ; 3 1 ' )  has  reported  a  large  variation  in  the 


solution  with  step -size,  as  can  be  seen  in  Fig  3.14,  the  deviation  has  not  been 
significant  in  the  results  obtained  using  .YOFEAP.  Fig 


snows  tr.e  piastic  Parceling  curve  :or  two  elements  in 
found  to  exactly  follow  the  specified  material  properties,  indicating  the  accuracy 
of  the  stress  reduction  scheme  in  NOFEAP.  The  bulging  behavior  of  the  cylinder 
has  been  compared  to  Yamaca's  solution  in  Fig  3.16,  and  the  agreement  is  found 
to  be  excellent. 


The  sample  examples  given  so  far  have  verified  the  elastic-plastic  and  geomet¬ 
rical!-/  non-linear  formulations  in  NOFEAP.  The  next  two  examples  compare  the 
contact  problem  formulation  with  the  Hertzian  theory. 


!.o  Hertz  Contact  Problem 


The  classical  Hertzian  contact  problem  is  selected  to  verify  the  accuracy  of  the 
contact  problem  formulation.  The  problem  involves  contact  between  two  cylinders, 
whose  longitudinai  axes  are  parallel.  The  two  cylinders  are  assumed  to  be  identical 
in  yeometrv  and  material  properties  in  the  following  examples.  This  assumption 
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computationally  simplifies  the  problem  by  symmetry  considerations  to  that  of  one 
cylinder  pressed  on  a  rigid  half-space. 

The  finite  element  mesh,  shown  in  Fig  3.17,  consists  of  117  elements,  each  with 
four  or  five  nodes.  The  details  of  region  ABCD,  near  the  contact  surface,  are 
given  in  Fig  3.18.  Only  a  portion  of  the  cylinder  has  been  modelled.  The  rigid 
half-plane  is  modelled  by  18  fixed  nodes.  Candidate  contact  nodes  are  distributed 
along  the  contact  surface  of  the  cylinder.  The  cylindrical  surface  has  a  radius  of 
curvature  of  254  mm.  A  single  concentrated  load,  P,  is  applied  at  the  node  at  left 
top  corner,  as  shown  in  Fig  3.17.  The  load  is  increased  from  an  initial  value  of 
1779.2  N  to  a  maximum  of  8S96  N  in  four  equal  increments. 

The  theoretical  solution  relating  the  load  to  contact  length  and  stress  distribu¬ 
tion  is  obtained  from  the  reference  :Sll.  The  numerical  resuls  are  compared  with 
the  Hertzian  solutions  in  Fig  3.19  through  3.22.  The  effect  on  contact  length  and 
maximum  nodal  stress  (traction)  at  various  loads  is  shown  in  Fig  3.19  and  3.20 
respectively.  The  numerical  results  differ  slightly  from  the  analytical  ones.  The 
distribution  of  stress,  however,  is  found  to  match  the  Hertzian  stress  distribution, 
as  seen  in  Fig  3.20. 

The  tangential  and  normal  stress  along  the  load  axis  are  plotted  in  Fig  3.21  and 
3.22  respectively.  The  stresses  drop  from  a  maximum  at  the  contact  surface  to  a 
negligible  value  at  a  depth  of  the  magnitude  of  the  contact  length.  These  stresses 
also  follow  the  Hertzian  distribution  but  differ  in  magnitude  by  a  small  factor. 
The  error  in  the  contact  length  and  stresses  is  found  to  be  less  than  5%  and  is 
considered  to  be  satisfactory. 
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3.6  An  Indentation  Problem 

The  indentation  of  a  cylindrical  punch  into  a  flat  semi-infinite  slab  provides 
a  variety  of  examples.  For  a  frictionless  interface  between  the  punch  and  an 
elastic  slab,  the  problem  can  be  solved  using  Hertz’s  contact  theory.  The  effect  of 
friction  along  the  interface  can  be  studied  for  indentation  into  an  elastic  or  elastic- 
plastic  slab.  The  problem  has  also  been  analyzed  in  detail  by  Cheng  [31].  Various 
combinations  of  material  and  interface  properties  have  been  used  to  provide  a 
number  of  solutions  to  the  indentation  problem.  Some  of  these  solutions  are 
presented  next. 

The  same  finite  element  mesh  (with  4-ncde  isoparametric  elements'!,  Fig  3.23. 
has  been  used  to  soive  the  different  cases  of  the  problem.  For  an  elastic  analysis, 
the  yieid  stress  is  set  to  an  appropriately  high  value.  The  siab  is  deformed  by 
pushing  the  cylindrical  punch  into  the  siab  in  steps  of  0.025  cm  until  a  maximum 
indentation  equal  to  half  the  thickness  of  the  siab  is  reached. 


The  elastic  deformations  and  the  distribution  of  effective  stresses  in  the  slab 
are  plotted  in  Fig  3.24  and  3.25  respective:;/.  Simiiar  plots  for  an  elastic-plastic 
analysis  with  a  non-zero  friction  coefficient  (0.3)  are  given  in  Fig  3.26  through  3.2S. 
The  main  difference  in  the  deformation  patterns  can  be  seen  along  the  far  edge  of 
the  siab.  While  the  elastic  siab  produces  more  lifting  behavior  and  less  expansion 
of  the  edge,  the  elastic-plastic  siab  produces  more  lateral  plastic  flow  than  the 
vertical  movement  along  the  edges.  That  there  is  significant  plastic  deformation 
in  the  latter  case  can  be  confirmed  by  observing  Fig  3.27,  whicn  suggests  that 
better  than  75?o  of  the  siab  has  permanently  deformed.  The  dirterence  in  tne 
effective  stress  contours  (Fig  3.25,  3.2$)  in  the  two  cases  can  a. so  be  attributed  t 
redistribution  of  stresses  due  to  plastic  fiow  in  the  eiastic-piastic  siab. 
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CONTACT  OF  RIGID  CYLINDER  WITH  PLASTIC  BLOCK 


ONjyj  i 


The  variation  in  maximum  contact  stress  (traction)  at  nodes  with  increasing 
depth  of  indentation  is  plotted  in  Fig  3.29.  The  stresses  have  been  normalized 
with  the  maximum  normal  stress  for  the  frictionless  case.  Although  the  difference 
in  frictional  and  frictionless  elastic  cases  is  not  significant,  all  of  these  produce 
considerably  higher  stresses  compared  to  the  elastic-plastic  slab.  The  change 
in  contact  zone  with  depth  of  indentation  is  shown  in  Fig  3.30.  The  effects  of 
friction  values,  indentation  depths  and  work-hardening  properties  of  the  material 
on  contact  conditions  have  been  plotted  in  Fig  3.31.  The  contact  stresses  have 
been  normalized  with  the  maximum  stress  for  the  first  curve  in  the  plots.  A  high 
slope  of  hardening  curve  produces  a  significant  increase  in  both  contact  stresses 
and  contact  length.  Increasing  the  friction  coefficient  along  the  interface  produces 
higher  stressess  without  affecting  the  contact  zone. 

The  contact  zone  for  plastic  contact  is,  expectedly,  higher  than  that  for  elastic 
contact.  The  frictioniess,  elastic  solutions  are  plotted  with  the  Hertzian  solutions 
[82]  in  Fig  3.32.  The  agreement  between  the  numerical  and  analytical  contact 
lengths  and  stresses  is  fairly  good.  The  solutions  of  Cheng  [31]  have  not  been 
reproduced  here  for  a  number  of  reasons.  Cheng’s  solutions  have  been  found  to 
differ  from  the  computed  NOFEAP  (and  Hertzian)  solutions  by  a  significant  factor. 
However,  Cheng  has  produced  plots  of  contact  stresses  versus  contact  length  which 
shows  an  exact  match  with  Hertzian  solution,  even  for  contact  lengths  of  the  order 
of  65$o  of  the  punch  radius.  According  to  the  hypothesis  of  Hertz’s  theory  [82!, 
the  theoretical  derivation  may  be  used  provided  that  the  surface  of  contact  is  very 
small  compared  to  the  radii  of  curvature  of  the  contacting  bodies  and  as  long 
as  the  surface  of  contact  is  along  a  plane  tangent  to  the  two  contacting  bodies. 
Clearly,  both  these  constraints  are  no  longer  satisfied  for  large  contact  zones  in 
the  indentation  problem.  With  these  observations  in  mind,  Cheng’s  solution  is 


not  consistent  with  the  Hertzian  theory,  and  therefore,  those  results  have  not 
been  included  for  comparison  in  this  thesis. 


The  examples  given  in  the  previous  sections  compare  the  non-linear  formula¬ 
tions,  individually  and  in  combinations,  with  available  numerical,  analytical  and 
experimental  calculations.  Several  other  problems,  such  as  stretching  of  a  necked 
or  perforated  specimen,  large  deflection  of  a  spherical  shell,  elastic  contact  between 
a  crankshaft  and  a  pin,  etc.,  have  been  solved  using  NOFEAP  and  the  results  have 
been  found  to  be  satisfactory.  The  next  example  combines  all  three  non-linearities 
into  one  problem  to  demonstrate  the  application  of  NOFEAP  to  a  metal  forming 
process. 

3.7  Head  Forming  Process 

The  head  forming  process  is  characterized  by  large  plastic  deformations  coupled 
with  contact  conditions  at  the  die-workpiece  interface.  The  finite  element  mesh, 
shown  in  Fig  3.33,  is  used  for  two  separate  cases  with  both  frictionless  and  perfectly 
sticking  friction  conditions  at  the  top  surface.  The  frictionless  case  also  assumed  a 
frictionless  die,  while  non-zero  friction  coefficient  has  been  used  at  the  horizontal 
die  with  no-slip  conditions  at  the  top  of  cylinder.  The  head  is  formed  by  reducing 
the  height  of  the  cylinder  in  steps  of  0.25%  of  the  initial  height  for  160  steps.  The 
material  properties  of  the  cylinder  are  assumed  to  be  the  following:  modulus  of 
elasticity  =  200  GPa,  yield  stress  =  700  MPa,  secondary  modulus  =  300  MPa, 
and  Poisson's  ratio  =  0.3. 

The  deformed  shape  of  the  cylinder  for  no-stick  and  no-slip  conditions  at  the 
top  interface  are  shown  in  Fig  3.34  and  3.35  respectively.  In  the  frictionless  case, 
the  top  surface  expands  laterally  indicating  that  the  punch  radius  is  assumed  to 
be  much  larger  than  the  radius  of  the  cylinder.  For  the  no-slip  case,  however, 


the  punch  and  the  cylinder  are  assumed  to  be  of  the  same  radius,  permitting  the 
bulging  of  the  cylinder,  as  shown  in  Fig  3.25. 

The  load-deflection  curve  for  the  two  cases  are  shown  in  Fig  3.36.  The  load  is 
computed  from  the  reactions  at  the  top  of  the  cylinder.  The  two  curves  separate 
when  the  cylinder  first  makes  contact  with  the  horizontal  part  of  the  die.  The 
frictionless  cylinder,  as  expected,  requires  less  force  to  reach  a  deformation  state 
than  the  rough  cylinder.  Finally,  the  contact  forces,  for  the  two  cases  at  the 
horizontal  part  of  the  die,  are  plotted  in  Fig  3.37.  The  friction  force,  at  the 
horizontal  die,  for  the  no-slip  case  causes  excess  force  to  be  required  as  shown  in 
Fig  3.36. 

The  examples  given  in  the  previous  sections  verify  various  non-linear  formula¬ 
tions  in  NOFEAP  and  also  demonstrate  the  applicability  of  the  program  to  metal 
forming  processes.  The  next  chapter  presents  a  model  of  the  rolling  process  and 
presents  the  analyses  of  deformations  in  rolling  cylinders  and  gear  teeth. 
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Chapter  4 

MODELLING  THE  AUSROLLING  PROCESS 

Previous  chapters  have  discussed  the  characteristics  of  ausrolling  and  the 
mathematics  of  the  non-linear  formulations  necessary  to  represent  the  plastic 
flow  in  ausrolling.  This  chapter  addresses  the  remaining  critical  feature,  from 
a  numerical  point  of  view,  of  ausrolling  which  is  the  rolling  process.  A  model  of 
the  rolling  process,  included  in  the  NOFEAP  program,  is  explained  in  the  next 
sections.  The  effectiveness  of  the  solution  algorithm  for  modelling  deformations  in 
a  disk  subjected  to  rolling  loads  is  compared  to  some  experimental  observations  in 
the  subsequent  sections.  The  deformations  in  a  gear  tooth,  when  the  gear  is  rolled 
against  a  hard  die,  are  modelled  using  NOFEAP.  These  results  are  presented  in 
the  last  section  of  this  chapter. 

4.1  A  Model  of  the  Rolling  Process 


A  numerical  model  of  the  sheet  metal  rolling  processes  has  been  attempted 
by  several  authors.  A  method  of  using  estimated  pressure  values  at  the  roll-sheet 
interface  [9]  has  been  popular  due  to  its  simplicity.  Another  approach  by  Bhargava 
[10,11]  simulates  rolling  by  translating  a  semi-elliptic  pressure  distribution  along 
the  surface  of  a  flat  workpiece.  The  advantage  offered  by  the  simplicity  of  the 
method  is  lost  when  an  attempt  is  made  to  appiy  the  method  to  modelling  plastic 
deformations  of  a  curved,  possibly  complex  surface.  Therefore,  the  rolling  process 
is  modelled  in  NOFEAP  using  some  geometric  and  kinematic  considerations. 

Consider  two  cylinders,  of  radii  Rq  and  Rq  respectively,  tangent  to  a  common 
plane  as  shown  in  Figure  4.1.  The  subscripts  D  and  G  refer  to  the  die  and  the 
deformable  body,  i.e.,  the  workpiece  gear  respectively.  The  two  circles  of  Figure  4.1 
may  be  taken  to  represent  the  pitch  circles  of  two  mating  gears,  and  the  cylinders 
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will  henceforth  be  referred  to  as  the  gear  and  the  die.  Assuming  no  slipping  along 
the  surface,  the  speeds  of  rotation,  fi d  and  tie,  of  the  die  and  the  gear  are  related 
by  the  following  constraint: 


Op  *  Rd  =  Hg  *  Rg‘ 


(4.1) 


so  that  the  distance  travelled  along  the  circumference  of  the  two  cylinders  is  the 
same  at  any  time.^  If  the  rotational  degree  of  freedom  (DOF)  of  the  gear  G  is 
constrained  while  the  translational  as  well  as  the  rotational  DOF  of  the  die  D  is 
unconstrained,  a  rotation  of  the  die  with  no  associated  slip  will  also  result  in  a 
translation  of  D  along  the  circumference  of  the  gear  resulting  in  a  planetary  gear 
motion  (Figure  4.2).  If  the  angle  travelled  along  the  circumference  of  the  gear  is 
taken  to  be  8q,  the  die  undergoes  two  mutually  separable  motions:  (l)  a  rigid 
body  translation  of  the  gear  as  it  is  rotated  about  the  center  of  the  gear  by  an 
angle  Qg »  and  (2)  a  rigid  body  rotation  of  the  die  about  its  own  center  by  an  angle 
6 d  given  by 

'Rd  +  Rg' 


eD  =  (— 


(4.2) 


This  concept  is  used  in  NOFEAP  to  model  the  rolling  process.  The  gear  is 
assumed  to  be  non-rotating  while  the  die  revolves  around  its  own  axis  as  well 
as  the  axis  of  the  gear.  The  model,  although  exceedingly  simple  in  concept,  has 
provided  satisfactory  results  with  the  total  angle  of  rotation  6g  applied  in  small 
enough  increments  as  shown  in  the  following  sections. 

Since  the  angle  of  secondary  rotation,  #£>,  is  a  function  of  the  radii  of  the  gear 
and  the  die,  by  changing  the  radii,  rolling  about  varying  surfaces  can  be  simulated. 


■ 


REPRESENTATION  OF  ROLLING  CYLINDERS  IN  “NOFEAP" 


Another  simple  modification  could  be  the  introduction  of  slipping  at  the  interface 
by  defining  the  secondary  angle  by 


9d  =  P(  -G)Qg  +  (1  -  (3)0S-  (4.3) 

where  0  <  fi  <  1  can  be  used  to  allow  conditions  ranging  from  pure  rolling 
(P  =  1)  to  pure  slipping  (/3  =  0),  and  9s,  the  angle  of  slipping,  taken  as  an 
independent  parameter.  This  correction  has  not  been  introduced  in  NOFEAP, 
since  pure  rolling  condition  at  some  radius  is  ensured  in  gear  rolling. 

This  model  of  the  rolling  process,  coupled  with  the  non-linear  formulations  of 
Chapter  2,  has  been  used  to  analyze  deformations  in  a  disk  subjected  to  rolling 
loads.  This  cumulative  model  is  compared  to  experimental  observations  as  given 
in  the  next  section. 

4.2  Study  of  Disk- Rolling 

In  order  to  test  the  model  of  rolling  process,  the  numerical  results  have  been 
compared  with  experimental  tests  of  rolled  aluminum  disks.  The  next  section 
describes  the  experimental  procedure  which  is  followed  by  a  description  of  the 
numerical  simulation  of  disk-rolling. 

4.2.1  Experimental  Procedure  and  Results 

A  number  of  experimental  rolling  tests  have  been  conducted  on  Aluminum  6061 
disks  using  a  round-faced  carbon  steel  die  and  a  flat-faced  high  speed  tool  steel  die 
(Figure  4.3).  The  inner  and  outer  radii  of  the  workpieces  were  0.701  and  1.7145 
inches  (17. SO  and  43.55  mm)  respectively,  and  their  thickness  was  one  inch  (25.4 
mm) .  The  outer  radius  of  the  round-faced  die  was  4.5905  inch  (116.60  mm)  and 
that  of  the  flat-faced  die  was  4.5S85  inch  (116.55  rum). 


The  ausrolling  machine,  whose  schematic  representation  is  shown  in  Figure 
1.3,  was  used  to  conduct  the  tests.  The  specimen  was  mounted  in  place  of 
the  gear  (B),  and  was  only  pushed  horizontally  into  the  die  (A)  by  moving  the 
rear  ram  downwards.  The  die  (A)  was  rotated  for  a  specified  time  as  soon  as 
the  workpiece  was  indented,  thus  working  a  portion  of  the  circumference  of  the 
latter.  The  parameters  controlling  the  depth  of  indentation  and  the  time  of 
working,  along  with  other  essential  parameters  were  set  and  controlled  by  an  8088 
microprocessor  based  personal  computer  and  interfaced  with  the  servo-controlled 
forming  equipment  using  a  Daytronic  System  10  data  acquisition  and  process 
control  system. 

The  primary  aim  of  the  experiments  was  to  compare  the  experimentally  and 
numerically  obtained  deformation  patterns.  The  force  values  obtained  from  the 
experiments  could  not  be  used  directly  as  the  force  was  the  resultant  force  required 
to  move  the  workpiece  and  to  keep  it  pressed  against  the  rotating  die.  Various 
deformation  patterns  were  obtained  by  indenting  the  workpiece  using  either  the 
round  or  the  flat  die  and  by  rotating  the  die  through  different  angles.  Three 
indentations  were  made  on  a  workpiece  using  the  round  die  while  only  one  test 
could  be  performed  on  a  piece  using  the  flat  die.  Various  indentation  depths  were 
explored  to  obtain  sufficient  or  measurable  residual  deformation  on  the  disks. 
The  complete  experiment  included  the  following  steps:  (l)  set  the  parameters  for 
indentation  depth  and  time  of  rolling;  (2)  feed  the  disk  into  the  die;  (3)  start 
rotating  the  die;  (4)  retract  the  disk  after  the  preset  time;  (5)  repeat  steps  1 
through  4  for  indenting  another  part  of  the  disk  (for  round  die  only);  and/or  (6) 
extract  the  disk  and  measure  the  residual  deformations. 

The  residual  deformations  were  measured  using  a  dial  guage  of  least  count 
0.0005  inch  (0.0127  mm).  One  undeformed  aluminum  disk  was  used  as  reference 


to  measure  both  the  initial  and  the  deformed  shapes  of  the  remaining  disks 
in  order  to  maintain  consistency  of  the  zero  setting  of  the  dial  guage.  The 
residual  deformations  (dotted  lines)  of  the  disks  were  plotted  by  superimposing 
the  magnified  deformations  on  the  initial  undeformed  shapes  (solid  lines)  of  the 
disks.  Some  of  these  plots  are  shown  in  Figures  4.4  through  4.10.  Also  displayed 
on  the  plots  is  the  information  such  as  the  die  type,  angle  of  rolling,  maximum 
force  required  in  the  experiment,  maximum  residual  deformation  as  indicated  by 
the  indentation  value,  etc.  The  induced  deformations  using  the  fiat  die  are  quite 
small  as  compared  to  those  produced  by  the  round  die  whereas  the  forces  follow 
an  opposite  trend.  Such  differences  in  magnitudes  is  to  be  expected  since  the  fiat 
die  works  the  total  thickness  of  the  disk  while  the  round  die  only  makes  a  thin 
indentation  of  the  order  of  ^  inch  (3.17  mm). 

Numerical  simulation  of  rolling  aluminum  disks  were  performed  using  NOFEAP 
to  reproduce  some  of  the  experimental  results.  A  discussion  of  the  development 
of  the  solution  algorithm  is  given  in  the  next  section. 

4.2.2  Numerical  Simulation 

Numerical  simulation  of  disk-rolling  consisted  of  one  initial  run  of  the  program 
to  roll  the  disk  by  one  complete  revolution.  Due  to  the  computation  time  required 
for  each  run,  it  was  not  possible  to  individually  simulate  each  experiment.  The 
force  and  the  deformations  were  obtained  at  a  number  of  intermediate  rolling  steps 
so  that  these  values  could  be  compared  to  some  experimental  results. 

The  input  data  file  used  in  the  analysis  is  given  in  the  Appendix  B.  Salient 
features  of  the  input  file  are  discussed  here.  Details  of  the  commands  used  can  be 
obtained  from  the  NOFEAP  user’s  manual  [81.  The  program  is  activated  by  the 
“XOFE”  macro  command  placed  at  the  top  of  the  file.  The  control  parameters 
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follow  the  NOFE  command.  The  first  line  indicates  the  number  of  the  bodies  in 
the  problem  and  is  set  to  two  in  order  to  accomodate  the  workpiece  (body  #1) 
and  the  die  (body  #2).  The  workpiece  is  modelled  by  144  nodes  and  96  elements, 
each  defined  by  five  or  less  nodes.  A  plane  strain  type  of  analysis  is  performed  and 
the  material  of  the  workpiece  is  assumed  to  follow  the  von  Mises  yield  criterion. 
The  geometric  and  material  description  follows  the  control  data.  The  properties 
of  Aluminum-6061  are  used  so  that  the  elastic  modulus  is  10,000  Ksi  (68.95  GPa), 
the  Poisson’s  ratio  is  0.33,  the  yield  stress  is  36  Ksi  (248.22  MPa)  and  the  slope 
of  hardening  curve  is  35.3  Ksi  (243,39  MPa).  Four  Gaussian  points  are  used 
in  each  element  for  integrating  the  stiffness  matrices  and  to  compute  the  strains. 
Geometric  non-linearities  are  not  included  in  the  analysis  to  keep  the  computation 
time  within  bounds.  The  workpiece  is  constrained  using  the  “BOUN”  command 
to  enable  it  to  translate  as  a  rigid  body.  The  die  is  only  represented  by  its  outer 
surface  and  is  modelled  by  129  nodes  listed  following  the  “BODY”  command. 

The  contact  surfaces  of  the  workpiece  and  the  die  are  defined  following  the 
“CONT”  command  in  the  mesh  input  phase.  The  contact  surface  of  the  former  is 
modelled  by  27  nodes,  which  represents  slightly  more  than  half  its  circumference, 
while  all  of  129  nodes  are  used  to  model  the  circumference  of  the  die.  Since  the 
number  of  degrees  of  freedom  in  the  contact  problem  solution  phase  is  a  function 
of  the  number  of  contact  nodes  on  the  workpiece,  but  not  of  those  on  a  rigid  die, 
including  large  number  of  nodes  on  the  latter’s  contact  surface  does  not  always 
significantly  increase  the  computation  time.  The  nodes  forming  the  contact  surface 
of  the  workpiece  are  listed  in  counterclockwise  direction.  The  coefficient  of  static 
friction  is  set  to  0.3.  The  solution  algorithm  follows  the  “MACR”  macro  command. 
The  comments  following  each  command  describe  the  function  of  the  command. 
Further  details  are  available  in  the  user’s  manual  [8] . 
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The  results  in  this  section  have  been  obtained  from  three  runs  of  NOFEAP  using 
three  meshes  of  the  workpiece  and  the  die.  The  first  run  rolled  the  workpiece  by 
360°  using  a  coarse  mesh  where  the  nodes  were  placed  at  an  angular  increment 
of  15°.  In  order  to  test  the  effect  of  the  mesh  size  on  the  solution,  a  finer  mesh 
was  placed  on  the  outer  surface  of  the  workpiece  with  an  angular  increment  of 
7^°  and  the  piece  rolled  by  180°.  Finally,  one-half  of  the  first  quadrant  of  the 
workpiece  was  again  subdivided  in  3-  increment  and  a  45  rolling  simulation  was 
performed.  The  results  of  these  numerical  simulations  are  given  in  Figures  4.11 
through  4.18.  Some  observations  on  the  similarities  and  the  differences  between 
the  numerical  and  experimental  results  are  included  next. 

4.3  Comparison  of  Numerical  and  Experimental  Results 

The  deformation  plots  at  intermediate  steps  and  the  final  deformed  shapes  of 
the  disks  for  the  experimental  observations  and  the  numerical  simulations  are 
compared  here.  It  should  be  noted  that  the  intermediate  numerical  plots  are  of 
deformations  while  the  die  is  still  in  the  process  of  rolling  the  workpiece  and 
hence  do  not  account  for  the  elastic  recovery  of  the  material.  On  the  other 
hand,  all  of  the  experimental  results  are  the  final  shapes  of  the  disks  after  the 
load  has  been  removed.  This  causes  certain  differences  in  the  magnitudes  of 
deformations.  However,  the  general  shapes  of  the  deformed  disks  compare  well 
for  the  experimental  and  the  simulated  disks  as  given  next. 

The  deformation  plots  are  grouped  according  to  their  angles  of  rolling  and  a 
comparison  is  attempted  as  follows: 

(l)  The  experimental  disk  3-3  (Figure  4.8)  with  rolling  angle  of  100°  can  be  matched 
against  the  numerical  plots  of  82p  and  127^°  (Figures  4.12,4.13).  The  fineness 
of  the  mesh  has  caused  the  sharp  change  in  the  deformation  at  the  end  of  rolled 
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zone  in  the  numerical  plots.  However,  the  tendency  of  the  material,  rolled  by 
the  die,  to  flow  over  the  original  shape  on  the  side  of  the  workpiece  opposite  to 
the  rolled  side  is  well  represented  in  the  numerical  plots. 

(2)  The  agreement  is  remarkable  for  the  experiment  number  1-2  (Figure  4.4)  and  the 
simulated  results  at  172^°  (Figure  4.14).  The  experimental  plot  is  considerably 
smoother  than  the  numerical  one  at  the  end  of  contact  zone  due  to  its  continuous 
nature  as  against  the  discreteness  of  the  finite  element  model.  The  material 
overflow  is  higher  for  the  numerical  model  but  this  could  be  due  to  some  lateral 
flow  in  the  experiments. 

(3)  The  rolled  disk  at  230°  (Figure  4.9)  is  stacked  against  Figure  4.15  (plot  at 
217g  ).  The  notable  exceptions  of  this  pair  from  the  earlier  ones  are  that  a 
coarser  mesh  is  employed  in  the  numerical  model  while  the  flat  die  is  used  in 
the  experiments.  Although  the  displacement  patterns  are  again  quite  similar, 
the  flat  die  does  not  produce  as  sharply  distinguishable  rolled  zone  as  brought 
forth  by  the  round  die.  The  reason  for  this  has  been  found  to  be  the  significant 
lateral  spread  of  the  material  in  the  former  case. 

(4)  The  plots  in  Figures  4.5  (experiment  2-1  at  2S0°)  and  4.16  (numerical  result  at 
262;;  )  compare  extremely  well  with  each  other  in  both  the  indented  and  the 
overflowing  material  zones.  The  resultant  cam  shaped  object  is  easily  identified 
in  the  two  plots. 

(5)  This  group,  formed  by  the  deformed  shapes  at  307 p  (Figure  4.17),  and 
experiment  2-2  (Figure  4.6)  with  the  round  die  and  number  6-1  (Figure  4.10) 
with  the  flat;  die,  both  with  rolling  angle  of  320°,  makes  some  interesting  points. 


The  numerical  model  again  predicts  a  cam  shaped  object,  as  also  produced 
by  the  round  die.  However,  the  flat  die  produces  very  insignificant  amount  of 
overflowing  material  while  causing  remarkable  indented  surface.  One  reason 
for  this  behavior  may  lie  in  the  deformation  measurement,  possibly  due  to  the 
eccentricity  of  the  mounted  disk.  Since  the  deformations  are  of  the  order  of  only 
0.0005”  (0.013  mm),  only  a  slight  eccentricity  would  affect  the  results.  This  lack 
of  discernible  rolled  zone  was,  however,  observed  in  all  other  specimens  rolled 
with  the  flat  die.  A  measurement  of  the  lateral  dimensions  of  the  specimen 
revealed  lateral  flow  of  the  order  of  the  depth  of  indentation  in  these  pieces.  This 
contradicts  the  assumption  of  plane  strain  deformation  in  the  specimens.  The 
round  die,  while  indenting  only  about  (3.17  mm)  wide  material,  may  have 
produced  plane  strain  deformation  more  closely  since  the  deforming  material 
was  constrained  by  surrounding  material  thus  forcing  the  flow  to  occur  in  the 
plane  of  rolling.  This  may  explain  the  semblance  of  the  numerical  results  with 
the  round  die  experiments  rather  than  the  flat  die  specimens. 

(6)  The  preceding  discussion  may  also  be  applied  to  find  a  similarity  between 
Figures  4.10  and  4.18. 

Although  the  displacements  of  the  numerical  solutions  have  compared  well  with 
the  experimental  results,  a  comparison  of  the  forces  could  not  be  made.  As  seen 
in  the  Figures  4.4  through  4.10,  the  forces  varied  from  about  48,000  to  90,000  Lbs 
(0.2  to  0.4  MN)  for  the  round  die  and  from  100,000  to  200,000  Lbs  (0.445  to  0.S9 
MN)  for  the  flat  die.  Different  force  readings  were  obtained  for  same  settings  on 
the  machine  thus  the  reliablity  of  these  values  is  questionable. 

The  variation  of  the  total  forces  with  rolling  angle  for  the  numerical  simulations 
is  plotted  in  Figure  4.19  for  various  meshes.  The  difference  between  the  15°  mesh 


and  the  7^°  mesh  is  significant  and  is  attributed  to  the  coarseness  of  the  mesh 
coupled  with  the  stiff  behavior  of  the  elements  in  shear  and  bending.  Refining  the 
mesh  to  3-  does  not  produce  further  change  in  the  force  values  indicating  the 
convergence  of  the  solution  with  the  mesh  size. 

The  model  of  the  rolling  process  in  NOFEAP  is  found  to  compare  satisfactorily 
with  the  experimental  runs  as  shown  in  this  section.  The  applicability  of  NOFEAP 
to  modelling  deformations  in  a  gear  tooth  when  the  tooth  is  subjected  to  rolling 
loads  is  demonstrated  in  the  next  section. 

4.4  Example  of  Deformations  in  Gear  Rolling 

A  gear  is  only  represented  here  by  half  of  its  one  tooth  (Figure  4.20).  The  gear 
is  rolled  against  a  hard  (rigid)  die  also  modelled  by  half  of  a  tooth.  Two  simulation 
excercises  have  been  performed  with  clockwise  and  counterclockwise  rotations  of 
the  gear  (opposite  for  the  die)  thus  modeiling  the  deformations  in  the  approaching 
and  the  trailing  sides  of  the  tooth  respectively. 

The  gear  tooth  is  modelled  by  167  nodes  with  139  elements  each  containing 
four  or  less  nodes.  The  geometric  constraints,  imposed  on  the  tooth  and  shown  in 
Figure  4.20,  are  slightly  inaccurate  as  they  do  not  model  the  bending  of  the  tooth. 
In  later  simulations,  some  of  the  constraints  on  the  tooth  were  released  to  observe 
the  effect  of  bending  on  deformations.  This  influence  was  found  to  be  significant. 
These  results,  however,  are  not  included  here. 

The  gear  and  the  die  dimensions  and  other  reievent  details  are  given  next. 
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The  surface  of  the  teeth  have  been  generated  using  the  definition  of  the  involute. 
For  modelling  the  trailing  side  of  the  tooth,  the  die  is  placed  with  its  tip  touching 
the  base  of  the  tooth  of  the  gear,  then  pushed  into  the  gear  and  rolled  in  5° 
increments.  The  resulting  deformations,  stresses  and  plastic  flow  patterns  are 
shown  in  Figures  4.21  through  4.24.  The  deformations  follow  the  expected  pattern 
shown  in  Figure  4.25  with  the  material  accumulating  towards  the  tip  while  forming 
a  slight  cavity  near  the  pitch  line  of  the  gear.  The  stress  contours  of  Figure  4.23 
indicate  the  concentration  of  high  compressive  stresses  near  the  position  of  the  die. 
The  plastic  flow,  shown  in  Figure  4.24,  shows  the  top  of  the  tooth  with  permanent 
deformation  due  to  restricting  the  bending  of  the  tooth. 

A  similar  analysis  is  made  to  mode!  the  approaching  side  of  the  tooth  with 
the  gear  and  the  die  placed  as  shown  in  Figure  4.26.  The  die  is  rolled  in  steps 
of  3°  in  the  clockwise  direction.  The  deformed  shapes,  shown  in  Figures  4.27 
and  4.30,  again  follow  the  expected  patterns  (Figure  4.23)  with  the  material  now 
accumulating  away  from  the  tip  and  at  the  pitch  line.  Stress  contours  of  Figures 
4.28  and  4.29  show  the  motion  of  high  compressive  stress  layer  as  the  die  roils  over 
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The  plastic  zone  again  dominates  the  top  of  the  tooth  (Figure  4.31)  as  in  the 
previous  case. 


Trailing 

Side 


Approach  • 
Side 


^  Pitch  Line 


Figure  4.23  EXPECTED  CE.~0P.WA7  1  CMS  IN  A  GEAR  TOOTH 


DUE  TO  ROLLING 


The  preceding  examples  demonstrate  the  successful  application  of  NOFEAP 
to  modelling  deformations  on  both  sides  of  a  gear  tooth  as  might  be  expected 
in  aus roiling.  The  rolling  process  itself  is  most  satisfactorily  modelled  in  the 
program  as  shown  in  this  chapter.  This  concludes  the  modeiling  of  essential 
features  of  ausroiling.  An  approximate  analysis  of  the  ausrolling  process  can  now 
be  reliably  performed  using  the  non-linear  unite  element  analysis  program.  Further 
improvements  can  be  mace  and  additional  capabilities  can  be  introduced  in  the 
program  in  future.  Some  of  these  erfcrts  along  with  the  concluding  remarks  are 
included  in  the  next  chanter. 
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Chapter  5 

CONCLUSIONS  AND  FUTURE  WORK 

The  main  objective  of  this  thesis  has  been  to  implement  the  essential  features 
of  the  ausrolling  process  in  a  computer  program  to  enable  approximate  simulation 
of  the  deformations  in  a  gear  following  ausrolling.  This  objective  has  been 
successfully  achieved  with  the  steps  listed  in  the  following  concluding  remarks. 

5.1  Conclusions 

1.  A  non-linear  finite  element  analysis  program,  NOFEAP,  has  been  developed 
with  capabilities  to  model  non-linear  behavior  due  to  elastic-plastic  flow 
with  work  hardening,  large  displacements  and  strains,  and  frictional  contact 
conditions. 

2.  Each  of  these  implemented  non-linear  formulations  has  been  satisfactorily 
verified  by  comparing  them  with  existing  analytical,  numerical  and  experimental 
results. 

3.  A  model  of  the  rolling  process  has  been  implemented  in  the  computer  program 
NOFEAP. 

4.  Several  rolling  experiments  have  been  conducted  with  aluminum  6061  disks  and 
a  round-faced  carbon  steel  die,  and  a  flat-faced  tool  steel  die.  The  experiments 
have  also  been  simulated  using  the  NOFEAP  program. 

5.  Deformation  patterns  obtained  experimentally  and  numerically  show  excellent 
similarities.  The  results  of  the  round-die  experiments  compare  more  agreeaoly 
with  the  simulated  results  than  do  the  results  of  the  flat-die  experiments.  This 
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has  been  attributed  to  the  lateral  flow  in  the  specimens  rolled  with  the  fiat  die  j 

which  violated  the  assumption,  of  plane  strain  deformations.  I 

i 

6.  The  NOFEAP  program  can  be  used  to  model  deformations  in  gear  rolling  as  J 

1 

demonstrated  by  results  of  simulated  runs  on  gear  teeth.  | 

1 

5.2  Suggestions  for  Further  Developments  ; 

As  mentioned  earlier,  only  some  of  the  features  of  ausrolling  have  been 
implemented  in  the  program  NOFEAP  in  this  research.  There  is  scope  for  further  : 

development  of  NOFEAP  in  order  to  obtain  a  more  accurate  representation  of  > 

ausrolling.  Some  of  these  directions  are  indicated  next.  \ 

1.  Modelling  of  temperature  distributions  in  the  gear  and  the  effect  of  thermal 
stresses  on  the  deformations. 


2.  Representation  of  the  time-temperature-transformation  curve  of  the  material  in 
the  simulation.  This  would  include  changes  in  the  material  properties  and  their 
gradients  as  the  deformations  progresses. 


3.  Extension  of  the  model  to  three  dimensions  to  simulate  the  actions  of  the 
swaging  and  the  shearing  operations. 


4.  Miscellaneous  improvements  in  the  computation  process  to  reduce  the  time 
taken  to  complete  a  simulation.  This  would  range  from  increasing  the 
efficiencies  of  some  algorithms  in  the  program  to  implementing  the  program  on  a 
supercomputer  enabling  effective  utilization  of  the  vector  and  parallel  processing 
techniques  on  these  machines. 


5.3  Limitations  of  the  Model 


The  numerical  formulations  implemented  in  NOFEAP  are  restrictive  to  some 

extent  as  pointed  out  next. 

1.  The  elastic-plastic  formulations  assumes  an  isotropically  expanding  yield  sur¬ 
face.  For  a  more  accurate  analysis,  kinematic  or  mixed  hardening  conditions 
may  be  introduced. 

2.  The  use  of  additive  elastic  and  plastic  strain  increments  may  be  objected  to 
by  those  using  the  multiplicative  decomposition  of  the  deformations  gradient. 
However,  the  superiority  of  either  theory  over  the  other  is  yet  to  be  proved. 
Hence  the  formulation  in  the  program  is  considered  to  be  satisfactory. 

3.  Since  the  contact  conditions  are  defined  at  discrete  nodes,  some  error  is  usually 
introduced  due  to  relative  lateral  sliding  of  the  nodal  pair.  Although  these  nodal 
pairs  may  be  updated  in  NOFEAP  at  each  step  by  connecting  the  nodes  closest 
to  one  another,  a  separation  equal  to  half  the  distance  between  adjacent  nodes 
may  still  exist.  Since  these  nodes  are  usually  placed  in  close  juxtaposition,  this 
error  has  not  been  found  to  be  significant. 
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Appendix  A 

LOAD  FACTORS  IN  CONTACT  PROBLEM  SOLUTION 


Assuming  the  notations  in  section  2.3.2,  the  load  factors  for  transitions  from 
one  contact  state  to  another  are  given  here.  If  the  contact  status  of  an  element 
does  not  change,  the  load  factor,  Ct ,  is  equal  to  one. 

1.  Open  to  Sticking  Contact: 

Initial  Gap  ( 8n ) 


Rel.  Normal  Disp.  Inc.  (/n) 

2.  Open  to  Sliding  Contact:  This  transition  is  only  permitted  if  the  coefficient  of 
friction  is  zero,  and  if  there  is  not  relative  tangential  sliding  between  the  bodies 
during  the  load  increment.  If  these  conditions  are  met,  the  load  factor  is  given 
by  (A-l). 

.  ^ 

3.  Sticking  to  Open  contact:  This  transition  requires  special  attention.  Both 


Okamoto  (6j  and  Salamon  [77j  have  defined  the  load  factor  as 

Rn 


a  =  a2  =  — 


ARn 


(A-  2) 


However,  for  initially  sticking  contact,  the  frictional  force  increment  is  indepen¬ 
dent  of  ihe  normal  increment,  so  that  even  though  using  (A-2)  for  load  factor 


assures  that  the  total  normal  force  given  by 


R  =  Rn  -f  OcARt 


is  zero,  the  total  frictional  force  which  is 


[A  ~  3) 


Rt  -  Rt  - f 


■a* 

on 

&& 


(a  -  *) 


may  not  be  zero  for  certain  problems,  resulting  in  a  non-zero  traction  at  a  free 
surface.  Such  a  situtation  has  been  encountered  in  this  research.  To  avoid  this 
problem,  further  check  is  essential  as  given  next. 


Si  =  Signum(Rt  +  ARt) 
aZ  Sl  (ARt  ~  SxflARn) 


(A -5) 

(A -6) 


Then  a.  =  Ck3,  if  is  smaller  than  a 2)  and  the  new  contact  status  is  set  to 
“slip.”  If  this  condition  is  not  satisfied,  an  error  message  is  issued  and  further 
computations  are  terminated. 

4.  Sticking  to  Sliding  Contact:  Letting  S  =  Signum(Rt  —  ARt), 


a  =  a4  =  s 


{fj.Rn  -  sRt) 
(ARt  —  SflARn) 


(A -7) 


All  the  cases  listed  by  Okamoto  [6]  and  Salamon  [77]  can  be  reduced  to  this 
simple  expression. 

5.  Sliding  to  Open  Contact:  Since  the  frictional  force  is  a  dependent  variable,  a 
tractionless  free  surface  is  obtained  after  releasing  the  contact.  The  load  factijfr. 
is  given  by  equation  (A-2). 

6.  Sliding  to  Sticking  Contact:  The  load  factor  is  set  to  zero  [77]. 
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Appendix  B 

INPUT  DATA  FILE  FOR  SIMULATION  OF  DISK-ROLLING 


The  following  data  file  has  been  used  to  simulate  the  aluminum  6061  disk- 
rolling  experiments  using  the  non-linear  finite  element  analysis  program  NOFEAP. 
The  non-essential  entries  (comments)  in  the  file  have  been  printed  in  lower  case 
characters  while  all  the  required  commands  are  in  upper  case  letters.  Detailed 


explanation  of  each  of  the  following  commands  can  be  obtained  from  the  NOFEAP 
users’  manual  [8], 


NOFE  ELASTIC-PLASTIC  ROLLING  OF  A  DISK\ 

2 

144  96  12  '  2  5  1  2  3  0  2  2 

129  0  0  0 

*********** 


npoi 

nele 

nvfi 

ntyp 

nnod 

nmat 

ncri 

ELEMent  listin 

cr 

o 

1 

1 

l 

2 

14 

13 

2 

1 

2 

3 

15 

14 

3 

1 

3 

4 

16 

15 

4 

1 

4 

5 

17 

16 

5 

1 

5 

6 

18 

17 

6 

1 

6 

7 

19 

18 

7 

1 

7 

8 

20 

19 

8 

1 

8 

9 

21 

20 

9 

1 

9 

10 

22 

21 

10 

1 

10 

11 

23 

22 

11 

1 

11 

12 

24 

23 

12 

1 

12 

1 

13 

24 

13 

1 

27 

25 

13 

14 

26 

14 

1 

29 

27 

14 

15 

28 

15 

1 

31 

29 

15 

16 

30 

16 

1 

33 

31 

16 

17 

32 

17 

1 

35 

33 

17 

18 

34 

18 

1 

37 

35 

18 

19 

36 

19 

1 

39 

37 

19 

20 

38 

20 

1 

41 

39 

20 

21 

40 

21 

1 

43 

41 

21 

22 

42 

ntil  ndof  ndim 


aufiti'vK: 


68 

1  106 

59 

83 

130 

69 

1  59 

107 

131 

83 

70 

1  107 

60 

84 

131 

71 

1  60 

108 

132 

84 

72 

1  108 

61 

85 

132 

73 

1  61 

109 

133 

85 

74 

1  109 

62 

86 

133 

75 

1  62 

110 

134 

86 

76 

1  110 

63 

87 

134 

77 

1  63 

111 

135 

87 

78 

1  111 

64 

88 

135 

79 

1  64 

112 

136 

88 

80 

1  112 

65 

89 

136 

81 

1  65 

113 

137 

89 

82 

1  113 

66 

90 

137 

83 

1  66 

114 

138 

90 

84 

1  114 

67 

91 

138 

85 

1  67 

115 

139 

91 

86 

1  115 

68 

92 

139 

87 

1  68 

116 

140 

92 

88 

1  116 

69 

93 

140 

89 

1  69 

117 

141 

93 

90 

1  117 

70 

94 

141 

91 

1  70 

118 

142 

94 

92 

1  118 

71 

95 

142 

93 

1  71 

119 

143 

95 

94 

1  119 

72 

96 

143 

95 

1  72 

120 

144 

96 

96 

1  120 

49 

»  73 

144 

* 

*  * 

* 

* 

* 

ele 

mat  nl 

n2 

!  n3 

n4 

)ORdinates  of  nodes 

1 

0.00000 

0.70100 

2 

0.35050 

0.60708 

3 

0.60708 

0.35050 

4 

0.70100 

0.00000 

5 

0.60708 

- 

0.35050 

6 

0.35050 

- 

0.60708 

7 

0.00000 

- 

0.70100 

8 

-0.35050 

- 

0.60708 

9 

-0.60708 

- 

0.35050 

10 

-0.70100 

0.00000 

11 

-0.60708 

0.35050 

12 

-0.35050 

0.60708 

13 

0.00000 

0.95000 

n6  n7  n8  n9 


I 


BS 


*<. 

/V 
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14 

0.47500 

0.82272 

15 

0.82272 

0.47500 

16 

0.95000 

0.00000 

17 

0.82272 

-0.47500 

18 

0.47500 

-0.82272 

19 

0.00000 

-0.95000 

20 

-0.47500 

-0.82272 

21 

-0.82272 

-0.47500 

22 

-0.95000 

0.00000 

23 

-0.82272 

0.47500 

24 

-0.45700 

0.82272 

25 

0.00000 

1.20000 

26 

0.31058 

1.15911 

27 

0.60000 

1.03923 

28 

0.84853 

0.84853 

29 

1.03923 

0.60000 

30 

1.15911 

0.31058 

31 

1.20000 

0.00000 

32 

1.15911 

-0.31058 

33 

1.03923 

-0.60000 

34 

0.84853 

-0.84853 

35 

0.60000 

-1.03923 

36 

0.31058 

-1.15911 

37 

0.00000 

-1.20000 

38 

-0.31058 

-1.15911 

39 

-0.60000 

-1.03923 

40 

-0.84853 

-0.84853 

41 

-1.03923 

-0.60000 

42 

-1.15911 

-0.31058 

43 

-1.20000 

0.00000 

44 

-1.15911 

0.31058 

45 

-1.03923 

0.60000 

46 

-0.84853 

0.84853 

47 

-0.60000 

1.03923 

48 

-0.31058 

1.15911 

49 

0.00000 

1.45000 

0.00 

50 

0.37529 

1.40059 

15.00 

51 

0.72500 

1.25574 

30.00 

52 

1.02530 

1.02530 

45.00 

53 

1.25574 

0.72500 

60.00 

54 

1.40059 

0.37529 

75.00 

55 

1.45000 

0.00000 

90.00 

56 

1.40059 

-0.37529 

105.00 

57 

1.25574 

-0.72500 

120.000 

58 

1.02530 

-1.02530 

135.000 

59 

0.72500 

-1.25574 

150.000 

60 

0.37529 

-1.40059 

165.000 

61 

0.00000 

-1.45000 

180.000 

62 

-0.37529 

-1.40059 

195.000 

63 

-0.72500 

-1.25574 

210.000 

64 

-1.02530 

-1.02530 

225.000 

65 

-1.25574 

-0.72500 

240.000 

66 

-1.40059 

-0.37529 

255.000 

67 

-1.45000 

0.00000 

270.000 

68 

-1.40059 

0.37529 

285.000 

69 

-1.25574 

0.72500 

300.000 

70 

-1.02530 

1.02530 

315.000 

71 

-0.72500 

1.25574 

330.000 

72 

-0.37529 

1.40059 

345.000 

73 

0.00000 

1.71450 

0.000 

74 

0.44375 

1.65608 

15.000 

75 

0.85725 

1.48480 

30.000 

76 

1.21233 

1.21233 

45.000 

77 

1.48480 

0.85725 

60.000 

78 

1.65608 

0.44375 

75.000 

79 

1.71450 

0.00000 

90.000 

80 

1.65608 

-0.44375 

105.000 

81 

1.48480 

-0.85725 

120.000 

82 

1.21233 

-1.21233 

135.000 

83 

0.85725 

-1.48480 

150.000 

84 

0.44375 

-1.65608 

165.000 

85 

0.00000 

-1.71450 

180.000 

86 

-0.44375 

-1.65608 

195.000 

87 

-0.85725 

-1.48480 

210.000 

88 

-1.21233 

-1.21233 

225.000 

89 

-1.48480 

-0.85725 

240.000 

90 

-1.65608 

-0.44375 

255.000 

91 

-1.71450 

0.00000 

270.000 

92 

-1.65608 

0.44375 

285.000 

93 

-1.48480 

0.85725 

300.000 

94 

-1.21233 

1.21233 

315.000 

95 

-0.85725 

1.48480 

330.000 

96 

-0.44375 

1.65608 

345.000 

97 

0.18926 

1.43760 

7.500 

98 

0.55489 

1.33963 

22.500 

99 

0.88270 

1.15036 

37.500 

100 

1.15036 

0.88270 

52.500 

101 

1.33963 

0.55489 

67.500 

102 

1.43760 

0.18926 

82.500 

103 

1.43760 

-0.18926 

97.500 

104 

1.33963 

-0.55489 

112.500 

105 

1.15036 

-0.88270 

127.500 

Mr 


s»: 


& 


my 
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106 

0.88270 

-1.15036 

142.500 

107 

0.55489 

-1.33963 

157.500 

108 

0.18926 

-1.43760 

172.500 

109 

-0.18926 

-1.43760 

187.500 

110 

-0.55489 

-1.33963 

202.500 

111 

-0.88270 

-1.15036 

217.500 

112 

-1.15036 

-0.88270 

232.500 

113 

-1.33963 

-0.55489 

247.500 

114 

-1.43760 

-0.18926 

262.500 

115 

-1.43760 

0.18926 

277.500 

116 

-1.33963 

0.55489 

292.500 

117 

-1.15036 

0.88270 

307.500 

118 

-0.88270 

1.15036 

322.500 

119 

-0.55489 

1.33963 

337.500 

120 

-0.18926 

1.43760 

352.500 

121 

0.22379 

1.69983 

7.500 

122 

0.65611 

1.58399 

22.500 

123 

1.04372 

1.36020 

37.500 

124 

1.36020 

1.04372 

52.500 

125 

1.58399 

0.65611 

67.500 

126 

1.69983 

0.22379 

82.500 

127 

1.69983 

-0.22379 

97.500 

128 

1.58399 

-0.65611 

112.500 

129 

1.36020 

-1.04372 

127.500 

130 

1.04372 

-1.36020 

142.500 

131 

0.65611 

-1.58399 

157.500 

132 

0.22379 

-1.69983 

172.500 

133 

-0.22379 

-1.69983 

187.500 

134 

-0.65611 

-1.58399 

202.500 

135 

-1.04372 

-1.36020 

217.500 

136 

-1.36020 

-1.04372 

232.500 

137 

-1.58399 

-0.65611 

247.500 

138 

-1.69983 

-0.22379 

262.500 

139 

-1.69983 

0.22379 

277.500 

140 

-1.58399 

0.65611 

292.500 

141 

-1.36020 

1.04372 

307.500 

142 

-1.04372 

1.36020 

322.500 

143 

-0.65611 

1.58399 

337.500 

144 

-0.22379 

1.69983 

352.500 

n 

4 - X  — + 

4 - y - 4 

comment 
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0.100E+08 

0.330E+00 

0.360E+05 

0.353E+05 

Li 

0.100E+01 

0.200E+01 

6. 

mnat* 

1 

e 

nu  thick 

dense  yield 

hard 

s 

el  typ 

g-  P • 

r. 

BOUN 

y 

1 

11 

0.00000 

0.01500 

2 

11 

0.00000 

0.01500 

6 

3 

ii 

0.00000 

0.01500 

? 

4 

11 

0.00000 

0.01500 

5 

11 

0.00000 

0.01500 

6 

11 

0.00000 

0.01500 

7 

11 

0.00000 

0.01500 

8 

11 

0.00000 

0.01500 

hj 

9 

11 

0.00000 

0.01500 

*r 

10 

11 

0.00000 

0.01500 

g 

11 

11 

0.00000 

0.01500 

fe 

12 

11 

0.00000 

0.01500 

Ti« 

A 

n 

xy 

x  -  dis 

y  -  dis 

i 

END 

a 

BODY=2 

V 

COORdinates  of  die 

1 

0.00000 

1.71450 

0.000 

A 

2 

0.22424 

1.71998 

2.800 

$ 

3 

0.44795 

1.73641 

5.600 

a 

4 

0.67059 

1.76375 

8.400 

5 

0.89163 

1.80193 

11.200 

V 

V 

6 

1.11054 

1.85086 

14.000 

$ 

7 

1.32680 

1.91042 

16.800 

Si 

8 

1.53989 

1.98049 

19.600 

9 

1.74930 

2.06087 

22.400 

5 

10 

1.95454 

2.15139 

25.200 

? 

11 

2.15511 

2.25183 

28.000 

12 

2.35053 

2.36194 

30.800 

if" 

13 

2.54034 

2.48148 

33.600 

14 

2.72409 

2.61014 

36.400 

15 

2.90133 

2.74762 

39.200 

£ 

16 

3.07164 

2.89359 

42.000 

A 

17 

3.23462 

3.04772 

44.800 

f. 

A 

18 

3.38988 

3.20961 

47.600 

19 

3.53704 

3.37891 

50.400 

r 

.* 

«• 

K 

% 

.s 

;*'.v 


•fi 


JO 

3.67576 

3.55518 

53.200 

>1 

3.80570 

3.73802 

56.000 

12 

3.92655 

3.92700 

58.800 

13 

4.03803 

4.12165 

61.600 

>4 

4.13986 

4.32151 

64.400 

!5 

4.23181 

4.52611 

67.200 

!6 

4.31366 

4.73496 

70.000 

17 

4.38521 

4.94755 

72.800 

!8 

4.44628 

5.16339 

75.600 

!9 

4.49674 

5.38195 

78.400 

10 

4.53646 

5.60272 

81.200 

:i 

4.56535 

5.82516 

84.000 

12 

4.58334 

6.04875 

86.800 

13 

4.59050 

6.30500 

90.000 

!4 

4.58502 

6.52924 

92.800 

15 

4.56859 

6.75295 

95.600 

;6 

4.54125 

6.97559 

98.400 

17 

4.50307 

7.19663 

101.200 

!8 

4.45414 

7.41554 

104.000 

i9 

4.39458 

7.63180 

106.800 

iO 

4.32451 

7.84489 

109.600 

.1 

4.24413 

8.05430 

112.400 

,2 

4.15361 

8.25954 

115.200 

13 

4.05317 

8.46011 

118.000 

:4 

3.94306 

8.65553 

120.800 

15 

3.82352 

8.84534 

123.600 

16 

3.69486 

9.02909 

126.400 

i7 

3.5573S 

9.20633 

129.200 

3.41141 

9.37664 

132.000 

19 

3.25728 

9.53962 

134.800 

iO 

3.09539 

9.69488 

137.600 

il 

2.92609 

9.84204 

140.400 

12 

2.74982 

9.98076 

143.200 

13 

2.56698 

10.11070 

146.000 

14 

2.37800 

10.23155 

148.800 

>5 

2.18335 

10.34303 

151.600 

16 

1.98349 

10.44486 

154.400 

’  »T 

>  t 

1.77889 

10.53681 

157.200 

18 

1.57004 

10.61866 

160.000 

19 

1.35745 

10.69021 

162. SCO 

10 

1.14161 

10.75128 

165.600 

il 

0.92305 

10.80174 

163.400 

12 

0.70228 

10.84146 

171.200 

13 

0.47984 

10.87035 

174.000 

14 

0.25625 

10.88834 

176.800 

15 

0.00000 

10.S9550 

180.000 

If: 

m 
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66 

-0.22424 

10.89002 

182.800 

67 

-0.44795 

10.87359 

185.600 

68 

-0.67059 

10.84625 

188.400 

69 

-0.89163 

10.80807 

191.200 

70 

-1.11054 

10.75914 

194.000 

71 

-1.32680 

10.69958 

196.800 

72 

-1.53989 

10.62951 

199.600 

73 

-1.74930 

10.54913 

202.400 

74 

-1.95454 

10.45861 

205.200 

75 

-2.15511 

10.35817 

208.000 

76 

-2.35053 

10.24806 

210.800 

77 

-2.54034 

10.12852 

213.600 

78 

-2.72409 

9.99986 

'  216.400 

79 

-2.90133 

9.86238 

219.200 

80 

-3.07164 

9.71641 

222.000 

81 

-3.23462 

9.56228 

224.800 

82 

-3.38988 

9.40039 

227.600 

83 

-3.53704 

9.23109 

230.400 

84 

-3.67576 

9.05482 

233.200 

85 

-3.80570 

8.87198 

236.000 

86 

-3.92655 

8.68300 

238.800 

87 

-4.03803 

8.48835 

241.600 

88 

-4.13986 

8.28849 

244.400 

89 

-4.23181 

8.08389 

247.200 

90 

-4.31366 

7.87504 

250.000 

91 

-4.38521 

7.66245 

252.800 

92 

-4.44628 

7.44661 

255.600 

93 

-4.49674 

7.22S05 

258.400 

94 

-4.53646 

7.00723 

261.200 

95 

-4.56535 

6.78484 

264.000 

96 

-4.58334 

6.56125 

266.800 

97 

-4.59050 

6.30500 

270.000 

98 

-4.58502 

6.08076 

272.800 

99 

-4.56859 

5.85705 

275.600 

100 

-4.54125 

5.63441 

278.400 

101 

-4.50307 

5.41337 

281.200 

102 

-4.45414 

5.19446 

284.000 

103 

-4.39458 

4.97820 

2S6.800 

104 

-4.32451 

4.76511 

289.600 

105 

-4.24413 

4.55570 

292.400 

106 

-4.15361 

4.35046 

295.200 

107 

-4.05317 

4.149S9 

298.000 

108 

-3.94306 

3.95447 

300.800 

109 

-3.82352 

3.76466 

303.600 

110 

-3.69486 

3.58091 

306.400 

111 

■3.00 ( oS 

3.40367 

309.200 
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H 

112 

-3.41141 

3.23336 

312.000 

B 

113 

-3.25728 

3.07038 

314.800 

■ 

114 

-3.09539 

2.91512 

317.600 

■ 

115 

-2.92609 

2.76796 

320.400 

JR 

116 

-2.74982 

2.62924 

323.200 

H 

117 

-2.56698 

2.49930 

326.000 

K 

118 

-2.37800 

2.37845 

328.800 

K 

119 

-2.18335 

2.26697 

331.600 

P 

120 

-1.98349 

2.16514 

334.400 

|| 

121 

-1.77889 

2.07319 

337.200 

nj 

122 

-1.57004 

1.99134 

340.000 

si 

123 

-1.35745 

1.91979 

342.800 

ii 

124 

-1.14161 

1.85872 

•  345.600 

1 

125 

-0.92305 

1.80826 

348.400 

126 

-0.70228 

1.76854 

351.200 

£ 

127 

-0.47984 

1.73965 

354.000 

K 

128 

-0.25625 

1.72166 

356.800 

129 

0.00000 

6.30500 

center 

I 

n 

x— > 

- y - * 
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END 

CONT 

27 
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number 

NODEs 


of  nodes  on  contact 
on  contact  surfaces 


surface  of  ‘a’  and  ‘b’ 
of  workpiece  (a),  and  die  (b) 


133 

85 

132 

84 

131 

83 

130 

82 

129 

,  81 

128 

80 

127 

79 

126 

78 

125 

77 

124 

76 

123 

75 

122 

74 

121 

73 

144 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 

121 

122 

123 

124 

125 

126 

127 
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* 

* 

* 

* 
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Miscellaneous  data 

NORM=0.2S0E+00 

TANG=0.250E+00 

FRIC=0. 300000 

SMAX=0.10E+14 

CHAR=0.2500 

END 

ANGLe  of  normals  on  contact  surface  of  disk 

WORK:  1,27,  .8250E+02,  1,  .7500E+01 

END 

PRIN 

END 


MACR 

DATA 

TRAN 

FILE 

STEF 

1.00 

FILE 

SOLV 

1.00 

FILE 

STAT 

1.00 

FILE 

NODE 

1.00 

TOL 

1.00 

1.00 

INCR 

2.00 

LOOP 

FEED 

2.00 

DATA 

INCR 

INCR 

0.00 

SETL 

ROLL 

1.00 

2.00 

LOOP 

ROLL 

25.0 

FILE 

STIF 

3.00 

SDSS 

1.00 

1.00 

LOOP 

ITER 

15.0 

FILE 

STIF 

3.00 

FILE 

SOLV 

3.00 

FILE 

STAT 

3.00 

STAC 

15.0 

25.0 

ELEM 

15.0 

25.0 

comm 

.... 

LOOP 

CONT 

10.0 

CONT 

USOL 

cove 

NEXT 

CONT 

BACS 

15.0 

25.0 

FORM 

CONV 

TIME 

start  macro  instructions 

open  stiffness  file 
open  solution  file 
open  static  condensation  file 
open  nodal  file 

tolerance  value  for  convergence=l% 
initial  feed=(2.00)  (0.0150) 
start  feed 

input  data  for  increment 
default  to  data  incr  value 
set  “roll”  to  1  when  “feed”  is  2 
start  rolling  -  one  revolution 
rewind  stiffness  file 
compute  stiffnesses 
start  iterations 
rewind  stiffness  file 
rewind  solution  file 
rewind  static  condensation  file 
reduce  dof  in  system  of  eqns 
compute  norm  &  tang  gaps 
and  generate  contact  elements 
start  contact  solution 
generate  contact  stiffness  &  loads 
unsymmetric  system  solution 
convergence  check  of  contact 

back  substitute  into  global  system 
compute  stresses  A  equiv.  forces 
plasticity  convergence  check 
update  displacements 
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NEXT 

ITER 

REAC 

CONT 

1.00 

output  contact  reactions 

NODL 

6.00 

output  nodal  values 

ROLL 

NEXT 

ROLL 

-15.0 

roll  by  -15.0  degrees 

NODL 

NEXT 

END 

FEED 

0.00 

output  final  nodal  values 

TRAN 

0.00 

0.0150 

INCR 

-1.00 

back  off  workpiece  by  1  unit 

ROLL 

RAD 

NODE 

1.71450 

4.59050  radii  of  workpiece  &  die 

129.0  center  of  die 

WORK  0.00 

END 

0.00  center  of  workpiece  (x-y  origin) 

INCR 

STOP 

+1.5 

back  workpiece  off  by  1.5 

/?Sf 


